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INTRODUCTION 


This  report  summarizes  work  done  under  Contract  F04611-86-C-0008 
during  the  period  1  March  1986  through  30  September  1988.  The 
principal  goals  were  to  develop  state-space  models  and  computational 
algorithms  for  control  of  beam  and  plate  type  structures,  and,  more 
generally,  to  increase  the  understanding  of  the  basic  problems 
associated  with  this  development.  The  state-space  approach  is  based 
on  a  distributed  parameter  model  of  the  structure  that  includes  the 
fundamental  partial  differential  equations  without  modal  truncation. 

The  approach  is  to  use  basic  physical  principles  to  write  down 
the  governing  partial  differential  equations,  construct  a  state-space 
model  from  the  governing  equations,  formulate  the  optimal  control 
problem  in  terms  of  the  state-space  model,  develop  a  convergent  ap¬ 
proximation  scheme  and  conduct  numerical  experiments  to  test  the  method. 
The  basic  view  is  that  it  is  best  to  avoid  the  approximation  step  until 
it  is  required  to  compute.  This  process  is  carried  out  for  various 
beam  models,  with  and  without  actuator  dynamics,  and  for  a  simple  ■ 
rectangular  clamped  plate.  The  basic  beam  systems  are  examined  and 
then  the  effect  of  delays  are  studied. 

In  addition  to  the  basic  program  we  conducted  a  survey  of  joint 
dynamic  models  and  attempted  to  evaluate  the  feasibility  of  constructing 
a  comprehensive  computer  code  that  would  be  general  enough  to  handle 
more  complex  structures  composed  of  basic  beam  and  plate  elements. 


1 


This  report  if  divided  into  nine  parts.  In  the  next  section  we 
present  a  general  discussion  of  the  types  of  control  problems  and 
structures  we  considered.  In  the  section  EQUATIONS  OF  MOTION  we 
summarize  the  basic  governing  differential  equations  for  the  various 
structures,  using  Euler-Bernoulli  and  Timoshenko  beam  theories  and 
simple  plate  theory,  and  include  actuator  dynamic  equations.  In  the 
STATE  SPACE  MODEL  section  we  construct  the  appropriate  state-space 
models  by  presenting  the  state  space,  the  dynamic  state  operator  and 
the  input  operator  for  each  problem.  We  conclude  this  section  by 
noting  the  general  framework  that  covers  all  the  problems  considered 
in  this  study.  THE  CONTROL  PROBLEM  section  contains  a  review  of  the 
basic  distributed  parameter  Linear-Quadratic  Regulator  (LQR)  problem. 
The  output  operators  are  detailed  for  each  of  the  structural  control 
problems.  The  specific  form  of  the  optimal  gain  operator  is  presented 
for  each  beam  and  plate  problem.  The  NUMERICAL  PROCEDURES  section 
begins  with  a  review  of  the  basic  approximation  theory  needed  to  solve 
the  LQR  problems.  Also  included  in  this  section  is  a  brief  description 
of  the  particular  schemes  used  in  each  problem  and  a  comparison  of 
open-loop  systems  for  Euler-Bernoulli  and  Timoshenko  theories.  A 
survey  of  joint  modeling  is  given  in  the  JOINT  MODELING  section. 
Numerical  results  for  the  control  problems  are  presented  in  the 
NUMERICAL  RESULTS  section  and  the  results  of  our  feasibility  study 
are  summarized  in  SOFTWARE  ISSUES.  The  last  section  is  devoted  to 
closing  remarks  and  suggestions  of  problem  areas  that  need  further 
study . 
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DESCRIPTION  OF  BEAM  AND  PLATE  PROBLEMS 


We  consider  a  variety  of  mathematical  problems  in  which  the 
fundamental  objective  is  to  quickly  damp  vibrations  in  various 
structures.  We  concentrated  our  effort  on  two  basic  structures;  a 
beam  with  a  tip-body  at  one  end  and  a  clamped  rectangular  plate. 

The  beam  structure  consists  of  a  beam  with  a  tip-body  at  one 
end.  The  other  end  of  the  beam  is  attached  to  a  hub  which  is  free  to 
rotate  about  a  fixed  axis.  We  study  only  motions  in  the  plane  normal 
to  the  hub's  rotation  axis.  Four  versions  of  the  problem  are  studied 
in  all.  Firstly,  we  employ  two  distinct  theories  to  model  beam  motions. 
The  simplest  and  most  common  is  the  Euler-Bernoulli  theory  which  does 
not  account  for  strain  energy  due  to  shear.  In  this  model  the  elastic 
energy  is  solely  due  to  normal  stretching  and  compression  of  the  beam 
"fibers"  as  the  beam  assumes  a  deformed  shape.  A  more  complete  theory, 
the  Timoshenko  model,  includes  strain  energy  due  to  shearing  and 
kinetic  energy  due  to  rotation  of  the  beam  face  elements. 

The  other  important  feature  we  study  is  the  model  for  the  torque 
actuator.  In  the  simplest  case  one  imagines  an  ideal  actuator  that 
can  instantly  generate  the  amount  of  torque  required  by  the  control 
law.  An  alternative  model  Includes  a  second-order  actuator  modelled 
with  £  delay  differential  equation.  If  the  delay  terms  are  set  to 
zero  one  has  a  usual  second-order  actuator  model;  however,  the  effect 
of  delays  (due  to  computational  requirements,  for  example)  are  shown 
to  be  important . 
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The  two  beam  models  (Euler-Bernoulli  and  Timoshenko)  are  combined 
with  the  two  actuator  models  (ideal  and  second-order  delayed)  to 
produce  four  distinct  beam  models  to  be  studied.  In  each  case  the 
objective  is  to  derive  a  control  law  in  state-feedback  form.  The  pur¬ 
pose  of  the  control  law  is  to  enhance  the  decay-rate  for  the  unwanted 
structural  vibrations. 

The  plate  problem  considered  concerns  transverse  motions  of  a 
rectangular  plate  which  is  clamped  along  its  edges.  The  "control"  is 
produced  by  modulating  the  amplitude  of  a  distributed  force  with  fixed 
spatial  variation.  That  is,  the  force  is  given  by 

f(t,x,y)  -  u(t)  ij>(x,y) 

where  <p(x,y)  is  given  and  u(t)  is  the  control.  As  before,  the  objec¬ 
tive  is  to  formulate  a  state-feedback  control  law  that  will 
substantially  improve  the  rate  at  which  vibrations  decay. 
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EQUATIONS  OF  MOTION 


In  this  section  we  present  the  equations  of  motion  for  the  five 
problems  considered  in  this  report.  In  particular,  we  summarize  the 
equations  for  the  beam  with  tip-body  and  the  clamped  plate  for  the 
following  cases: 

•  An  Euler-Bernoulli  beam 

•  A  Timoshenko  beam 

.  An  Euler-Bernoulli  beam  with  actuator  dynamics 

•  A  Timoshenko  beam  with  actuator  dynamics 

•  An  isotropic  rectangular  plate 

A  detailed  derivation  of  the  Euler-Bernoulli  model  is  given  in  Ref.  2. 
Inclusion  of  actuator  dynamics  are  discussed  in  Ref.  2  and  the  extension 
to  the  Timoshenko  beam  model  follows  directly  from  Ref.  1,  Ref.  2  and 
Ref.  3.  The  reader  is  referred  to  Ref.  4  and  Ref.  5  for  detailed 
discussions  of  plate  modeling. 

BEAM  EQUATIONS 

We  present  the  various  mathematical  models  of  the  beam  structure 
shown  in  Figure  1.  It  is  assumed  that  the  structure  is  pivoting  about 
a  fixed  pivot  at  point  0,  and  that  the  motion  is  in  a  plane.  The 
development  is  achieved  by  direct  application  of  Newton's  Laws  to 
individual  members  of  the  structure  with  the  resulting  equations  summed 
to  determine  the  overall  motion. 

The  structure  consists  of  three  parts,  the  main  frame,  or  mass  A, 
the  beam,  or  mass  B,  and  the  end  mass,  or  mass  C.  The  main  frame  is 
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assumed  to  be  a  rigid  body  pivoting  about  point  0.  The  beam  is  assumed 
to  be  flexible  and  rigidly  attached  to  the  main  frame  as  a  cantilevered 
beam.  The  end  mass  is  rigidly  attached  to  the  end  of  the  beam  so  that 
it  moves  and  rotates  with  the  end  of  the  beam.  Finally,  the  attached 
point  is  assumed  to  be  at  the  center  of  mass  of  the  end  mass . 

The  coordinate  system  to  be  used  is  fixed  in  the  main  frame  and 
rotates  with  it.  The  origin  is  at  point  0  with  the  x  axis  pointing 
along  the  undeflected  beam,  the  z  axis,  the  axis  of  rotation,  and  the  y 
axis  completing  the  right  hand  set.  The  governing  equations  depend  of 
course  on  the  assumptions  about  the  beam.  We  summarize  these  equations 
for  the  four  cases  listed  above. 

Euler-Bernoulli  Beam  Theory 

It  follows  from  Ref.  2  that  the  motion  of  the  structure  is 


governed  by  the  system  of  equations 

0(t)  -  U)(t)  (1) 

I^,L(t)  -  El  u^(t,0)  +  M^(t)  (2) 

m_,[L(L(t)  +  n(t)l  -  El  u^^(t,L)  ,  (3) 

c  xxx 

I^[(L(t)  +  5(t)]  -  -El  u^^(t,L)  (4) 

u  (t,x)  +  Xw(t)  ■  -  —  u  (t,x)  (5) 

tt  p  xxxx  ' 


where  u(t,x)  denotes  the  deflection  of  the  beam,  e(t)  Is  the  angle  of 

hub  rotation,  n(t)  Is  the  relative  lineal  velocity  of  the  tip  mass, 

^(t)  Is  the  relative  angular  velocity  of  the  tip  mass  due  to  deflection 

and  M.(t)  Is  the  applied  moment  at  the  hub. 

A 
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Since  the  beam  remains  joined  to  the  hub,  one  has  the  canti¬ 
levered  boundary  conditions 

u(t,0)  ■  u^(t,0)  -  0 
u  (t,0)  -  u  (t,0)  «  0 

X 

and  the  integrity  of  the  upper  joint  requires  the  boundary  conditions 

u^Ct.L)  -  n(t) 

It  is  important  to  note  that  both  (6)  and  (7)  are  geometric  (i.e. 
essential)  boundary  conditions. 


Timoshenko  Beam  Theor\ 


The  Timoshenko  beam  model  includes  the  rotary  inertia  and  shearing 
deformations.  If  u(t,x)  denotes  the  deflection  of  the  beam,  then  the 
total  slope  u^(t,x)  is  given  by 


u^(t.x)  -  i|/(t,x)  +  6(t,x)  (8) 

where  'P(t,x)  is  the  slope  of  the  deflection  curve  when  shearing  is 
neglected  and  6(t,x)  is  the  angle  of  shear  at  the  neutral  axis.  The 
bending  moment  M(t,x)  and  shear  force  V(t,x)  are  given  by 


M(t,x)  «  El  '(’^(t,x) 

and 


(9) 


V(t,x)  -  -  K'AG  8(t,x)  -  -  K'AG  [u^(t,x)  -  ip(t.x)]  . 

respectively.  Here  A  is  the  cross-section  area,  K'  is  a  constant 
depending  on  the  shape  of  cross-section  and  G  is  the  modulus  of 
elasticity  in  shear  (Ref.  3). 
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The  equations  for  the  Timoshenko  beam  become 


(’(t)  =  w(t)  (11) 

I^;(t)  =  El  -y^Ct.O)  +  M^(t)  (12) 

m^[L,L(t)  +  n(t)]  =  -  K'AG  [u^(t,L)  -  ^(t.L)]  (13) 

I^[(L(t)  +  ^(t)]  =  -  El  ^^(t.L)  (14) 

(Utt(t,x)  +  X  w(t)]  =  K'G  [u^^(t.x)  -  i;^^(t,x)]  (15) 

pi  w(t)]  =  El  +  K'AG  [u^(t,x)  -  iKt,x)]  (16) 


Again,  one  has  the  cantilevered  boundary  conditions  at  the  hub,  x  =  0, 

u(t,0)  =  u  ^t,0)  =  0  1 

\  (17) 

Kt,0)  =  .i<  (t,0)  =0  J  . 


Also,  the  upper  joint  requires  the  boundary  conditions 

Uj.(t,L)  =  n(t) 

(;,^(t,L)  =  ^;(t) 


(18) 


Euler-Bernoulli  Beam  with  Actuator  Dynamics 

Returning  to  the  Euler-Bernoulli  model  (1)  -  (5)  we  add  a  system 
of  equations  that  incorporates  actuator  dynamics  with  a  delay.  Thus, 
the  total  system  for  the  Euler-Bernoulli  beam  with  tip-body  and 
delayed  actuator  dynamics  becomes 


i3(t)  =  a<t) 


(19) 


w(t) 


[EI^ 

u  (t,0)  + 

fl  ] 

.^A. 

XX 

.^A. 

[Cj  Xj(t)  +  Cn  X2(t)l 


(20) 
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/  ' 

[Lu(t)  +  n(t)]  -  ^  (21) 

in  XXX 
^  c 

[u(t)  +  l(t)]  -  [■  1^1  u^(t.L)  (22) 

[u  (t,x))  -  (u  (t,x)  +  xuj(t)]  (23) 

dt  XX  g^2  t 

^  [u  (t.x)  +  xu(t)]  -  [u  (t,x)]  (24) 

t  I  PJ  3x2  XX 

with  general  second-order  actuator  model 

Xl(t)  -  X2(t)  (25) 

X2(t)  -  aiixi(t)  +  ai2X2(t)  +  a2iXi(t-r)  +  a22X2(t-r) 

+  biu  (t)  ,  (26) 

^  c 


where  u  (t)  is  the  commanded  input  to  the  actuator  and  M  (t)  * 
c  A 

CjX^(t)  +  C2X2(t)  is  the  output.  The  time  delay  r  >  0  was  Included  to 
Illustrate  the  problems  one  can  encounter  if  numerical  schemes  for 
control  designs  are  improperly  constructed.  Note  that  if  a2i  ■  a22  “  0» 
then  (25)  -  (26)  is  a  second  order  linear  actuator. 


Timoshenko  Beam  with  Actuator  Dynamics 


If  one  rewrites  (11)  -  (16)  and  adds  (25)  -  (26)  to  these  equa¬ 


tions,  we  obtain  the  following  equations 


9(t)  -  w(t) 


El)  r 1 1 

a)(t)  -  —  tp  (t,0)  +  —  tcixi(t)  +  C2X2(t)] 

*  (V 

[Lw(t)  +  n(t)]  -  f'^^1  [u^(t,L)  -  iJ/(t,L)] 

la  X 

\  c  ' 
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[a,(t)  +  n(t)]  *  ^ 


-El 

I 


(30) 


lu^|.(t,x)  +  X  w(t)]  p 


K’G 


(31) 


['i'^j.(t,x)  +  0j(t)  ]  = 


E 

IPJ 


K  ’  Ap'^ 

-^yj  tu^(t,x)  -  ij>(t,x)]  (32) 


xi (t)  =  X2(t)  (33) 

X2(t)  =  aiixi(t)  +  ai2X2(t)  +  a2ixi(t-r)  +  a22X2(t-r) 

+  bi  u^(t)  ,  (34) 

where  again  M  (t)  =  cixi(t)  +  C2X2(t)  is  the  output  to  the  actuator 

A 


(33)  -  (34). 


Isotropic  Rectangular  Plate 

Figure  2  shows  a  plate  with  uniform  thickness  h  that  is  assumed 
small  in  comparison  with  its  other  dimension.  The  x-y  plane  is  taken 
to  be  the  (undeflected)  middle  plane  of  the  plate.  Deflections  in  the 
z  direction  are  assumed  small  in  comparison  to  the  thickness  and  the 
normals  to  the  x-y  plane  are  assumed  to  remain  normal  to  the  deflected 
middle  surface  during  vibrations.  The  governing  equation  for  the 
loaded  undamped  plate  becomes 


in  u^^  (t,x,y)  +  D42u(t,x,y)  =  f(t,x,y) 


(35) 


(  Eh 


,  V  is  Poisson's 


where  m  is  the  mass  per  unit  area,  D  * 
ratio  and  f(t,x,y)  is  an  applied  load.  Here  a  is  the  Laplacian  operator 


9x2  9y2 


so  that 


(36) 
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12 


3^*  2  3^  3** 

A2u(t,x,y)  =  —  u(t,x,y)  + — -  u(t,x,y)  +  -  u(t,x,y)]  .  (37) 

dx*  3x^9y^  Jy’* 

Let  U  =  (0,a)  ^  (0,b)  denote  the  surface  of  the  plate  and  ail 
denote  the  boundary  (i.e.  edges)  of  the  plate.  Then,  (35)  must  hold 
for  all  t  ^  0  and  (x,y)  6  and  the  clamped  boundary  conditions  are 
given  by 

u(t,x,y)  =  Au(t,x,y)  =  0,  (x,y)  e  (38) 

We  shall  consider  a  control  problem  for  the  plate  equation  (35)  with 
boundary  conditions  (38)  . 


13/14 


STATE  SPACE  MODELS 


In  this  section  we  construct  the  state-space  models  for  each  of 
the  problems  discussed  above.  These  models  are  used  later  to  formulate 
control  problems,  provide  explicit  representations  of  optimal  control 
laws  and  as  a  basis  for  constructing  computational  algorithms. 


EULER-BERNOULLI  BEAM 

In  order  to  construct  a  state-space  model,  equation  (5)  Is  first 
written  in  the  form 


±  tv''-**' 


(39) 


As  noted  In  Ref.  5,  this  system  suggests  that  the  functions  x  u  (t,x) 

XX 

and  X  ->■  [u^(t,x)  +  xu(t)]  are  natural  state  components.  Thus,  the 
"state  vector"  y(t)  is  defined  to  be 


y(t) 


9(t) 

U)(t) 

n(t)  +  La)(t) 
^(t)  +  (i)(t) 


u 

XX 


(t,x) 


U^(t,x)  +  XU)(t) 


The  appropriate  state  space  Is  the  Hilbert  space 


y  -  R  X  8  X  a  X  a  X  L2(0,L)  L2(0,L) 


(AO) 


(Al) 
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with  inner  product 


«y.y>Y  =  yiyi  +  y^y^  +  y  jVj  +  y^y, 


L  ^  .L 

El  y5(x)y5(x)dx  +  p  y6(x)yg(x)dx  . 

0  JO 


(42) 


The  space  'J  with  inner  product  (42)  is  a  Hilbert  space  and  ||y||  = 


/<y »y>jj  defines  a  norm  on  jj.  The  product  <y,y>y  is  essentially  the 
mechanical  energy  in  the  physical  system  at  state  y. 

Let  (Sp  denote  the  operator  that  acts  on  functions  and  produces 
the  value  of  that  function  at  point  p,  i.e.  6piJ>(*)  “  ij)(p)  •  Moreover, 

D  denotes  the  partial  differential  operator  .  Define  the  operator 


oX 


Aj.  on  D(i^)  c  y  by 


0  0 


0 

0 


0 

0 


0  0  0  0 


0  0  0  0 


0  0  0  0 


0  0  0  0  —  6,oD 

m  Lx 


n 

c 

-El 


T  "l 

c 

0 


D  2 

0  X 


0 

0 


(43) 


where  the  domain  of  A^,  is  given  by 

(yi .y2.y3.yi,.y5(*) .y6(*))^  e  U  /  y^.y^  h-(o,l), 
.  yo<o)  =  0  •  y6<o>  =  yi  •  =  y3»yb(^)  * 


D(Ag) 


(44) 


Let  B-,  be  the  input  operator  R -*•  U  defined  by 
E  E 


16 


0 

'a 
0 

0 

0 

0 

and  consider  the  system  in  y  given  by 

y(t)  -  Agy(t)  +  Wit)  (46) 

It  was  shown  in  Ref.  1  that  Ag  generates  a  CQ-semlgroup  on  y  and 

that  the  variation  of  parameters  formula 

y(t)  -  e^E*^y(0)  +  B  M  (s)ds  '  (47) 

Jq  ^  a 

provides  weak  solutions  of  the  partial  differential  equations  (1)  - 
(6)  with  initial  data  in  y.  In  particular  one  has  the  following  result. 

THEOREM  1.  If  y(0)  -ye  D(Ag)  and  M^(t)  belongs  to  Li(0,T), 
then  y(t)  defined  by  (47)  is  a  solution  to  the 
system  (1)  -  (5)  with  initial  data  yo* 

Theorem  1  provides  the  basic  theoretical  result  that  allows  one 
to  construct  convergent  numerical  schemes  for  control  design  (see  Refs. 
1,  2  and  6-15).  Later  we  augment  the  state  space  model  (46)  by 
adding  the  actuator  dynamics  with  delays. 
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It  is  important  to  understand  that  the  statement  that  yg  6  D(Ao) 
is  initial  data  for  (1)  -  (5)  means  that 

yg  =■  col  (0o,<iJo,ao»Yo*  VgCx)) 

satisfies  (44)  and  the  solution  of  (1)  -  (6)  satisfies  the  initial 
condition  0(0)  -  0o.  (^(0)  -  a)o,  [Lw(0)  +1(0)]  -  oq,  [u)(0)  +  C(0)]  »  yq. 
u^(O.x)  ■  Uq(x)  and  [u^(0,x)  +  xw(0)]  »  Vo(x). 

TIMOSHENKO  BEAM 

We  first  rewrite  the  second  order  equations  (15)  and  (16)  as  the 
first  order  system 

[u^(t,x)  +  xa)(t)]  -  ^  [u^(t,x)  -  <j^(t,x)]  (48) 

^  [*.^(c.x)  +«(t))  .  [ij  +  (£^}  [u__(t.x) 

-  •J»(t,x)]  (49) 

^  i|<^(t,x)  -  [i|/^(t,x)  +  a)(t)]  (50) 

'it  ~  '|j(t,x)]  -  [u^(t,x)  +  xa)(t)]  -  [i|;^(t,x)  +  u)(t)]  .  (51) 

Let  w(t)  denote  the  "state-vector" 
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w(t) 


0(t) 

l^(t) 

Lw(t)  +  n(t) 
w(t)  +  -;(t) 
u^(t»x)  +  xw(t) 
i/;^(t,x)  +  w(t) 

u^(t,x)  -  ii'(t,x) 


(52) 


The  appropriate  state  space  for  Timoshenko  beam  is  the  Hilbert  space 

U)  =  R'<RxRxRxl^xl2'<L2’'L2,  (53) 

where  L  2  L2(0,L),  and  the  inner  product  is  given  by 

<w,w>jjj  =  +  I^W2W2  “(,''3^3  +  I^W4W4 


+ 


L 

0 

•L 

•'0 


{Ap  W5(x)  WjCx)  +  Ip  Wg(x)  WgCx)}  dx 


{El  W7(x)  WyCx)  +  (K'AG)  Wg(x)  Wg(x)  }dx 


(54) 


Here  again  (54)  defines  an  inner  product  on  (53)  and  <w,w''|jj  =  llwll^ 
is  essentially  the  mechanical  energy  in  the  physical  system  at  state  w. 
Define  the  differential  operator  with  domain  D(iVp)  c  U)  by 
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0  10  0  0 


0  10  0  0 


0  0  0  0  0 


0  0  0  0  0 


0  0  0  0  0 


0  0  0  0  0  0 


0  0  0  0  0  D 


1_0  0  0  0  D  -1 


where  the  domain  of  is  given  by 


f-K'AGl 


K'AG' 


D(Aj.) 


J-  W  =  (WpW  ,,W3,W,,,Wr,(  •)  .Wg(  •)  »W7(  •)  »W2(  •) )  6  UJ/  I 

j  w.^(*).  Wg(*).  W7(*).  Wg(0  e  H^O.L)  and  ^ 

[  Wg(0)  =  0,  WgCD  =  Wj,  Wg(0)=W^,  ^  ' 


Let  be  the  input  operator  Rj,  !  R  -*•  IB  defined  by 

r »  1 
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and  consider  the  system  in  U)  given  by 


w(t)  =  Aj,  w(t)  +  Rj,  M^(0  . 


(58) 


The  following  result  may  be  established  using  the  methods  in  Ref.  1. 


THEOREM 


2.  The  operator  is  bounded  and  Aj.  generates  a  C^-semigroup 
e^*"  on  the  Hilbert  space  U). 


It  also  follows  that  if  Wg  €  D(Aj,)  then 


w(t)  =  Wq  +  M  (s) 

Jo 


(59) 


provides  weak  solutions  to  the  Timoshenko  beam  equations  (11)  -  (16) 
with  initial  data  Wg . 

We  turn  now  to  the  addition  of  the  actuator  dynamics.  The 
general  second-order  model  (25)  -  (26)  will  be  written  as  a  system  and 
added  to  the  state-space  models  for  the  Euler-Bernoulli  beam  (A6)  and 
the  Timoshenko  beam  (58),  respectively. 


EULER-BERNOULLI  BEAM  WITH  ACTUATOR  DYNAMICS 

T  XT 

Let  X  *  [xi,X2]^,  Cj  =  [Ci,C2]  ,  Bj  =  [O.b^]  and 


1 

0 

1  1 

o 

1  ^ 

o 

> 

11 

II 

< 

.  ^11 

^12  . 

1  1 

.  ®21 

®22  . 

The  system  (25)  -  (26)  can  be  written  in  matrix  form  as 
x(t)  =  Aix(t)  +  A^xft-r)  +  Biu^(t) 
and  the  applied  torque  M^(t)  is  given  by 

M^(t)  =  Cix(t) 


(60) 


(61) 
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The  addition  of  actuator  dynamics  to  the  state-space  model  for  the 
Euler-Bernoulli  beam  is  accomplished  by  combining  (46) ,  (60)  and  (61) 
into  "state-space  form" 

y(t)  -  Agy(t)  +  Bg  Cix(t) 
x(t)  -  Aix(t)  +  A2x(t-r)  +  Biu^(t) 

In  order  to  show  that  the  augmented  system  (62)  is  well-posed  we  must 
Introduce  an  augmented  state-space.  We  consider  two  distinct  cases, 
A2  *  0  and  A2  0. 


(63) 


(64) 


(65) 


(66) 


(67) 
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The  augmented  system  (62)  becomes 

zO(t)  =  Aq  z(t)  +  Bou^(t) 

Moreover!  one  can  easily  establish  the  following  result. 

THEOREM  3.  The  operator  Bg  is  bounded  and  Aq  generates  a  ( 
on  the  Hilbert  space  2°. 

The  case  with  delays  is  more  complex. 

Case  2:  Delayed  Actuator 

Let  Zg  denote  the  Hilbert  space 

X  L2(-r,0;n2) 

with  inner  product 

fO 

=  <z'^,z®>0g  +  <4i(s)  ,  (p(s)>  ds 

■'  -r 

‘T  ^  ^  iji 

where  z  =  (z'^.t))  €  and  z  =  (z°,4.)  6  respectively, 

on  D(A)  c  2g  by 

and 


r  1 

y 

■  V  +  BgCiX 

A 

X 

- 

A^x  +  A2<}>(-r) 

.  . 

4>'(-) 

The  operator  B  :  R  -►  2_  is  defined  by 

L 


D(A)  = 


y 

X 


(68) 

^-semigroup 

(69) 

(70) 

Define  A 

(71) 

(72) 
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0 


B  - 


0 


(73) 


The  following  result  is  non-trivlal.  However,  the  proof  can  be  given 
based  on  the  results  in  Ref.  1  and  Refs.  6-10. 


THEOREM  4.  The  operator  B  is  bounded  and  A  generates  a  C  -semigroup 

o 

At 

e  on  the  Hilbert  space  2^. 

It  follows  that  a  state-space  model  for  the  Euler-Bernoulli  beam  with 
delayed  actuator  dynamics  is  given  by 

i(t)  -  Az(t)  +  Bu^(t)  (74) 

c 

where  A  and  Bare  defined  by  (71)  -  (73).  Moreover,  for  initial  data 
Zq  *  (yo>*0»*^o)  ®  B(A)  and  u^(’)  6  L2(0,T),  (74)  has  a  unique  mild 
solution 

2(t)  -  e*^  Z0+  f  Bu  (s)  ds  (75) 

that  provides  a  weak  solution  to  the  partial-functional  differential 
equations  (19)  -  (26) . 


TIMOSHENKO  BEAM  WITH  ACTUATOR  DYNAMICS 

We  proceed  as  in  the  case  for  the  Euler-Bernoulli  beam.  If 
A2  “  0,  then  set 

(76) 

with  inner  product 
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(77) 


<Zo*z°>Ot  “  ’'2^2 


Let  D(An)  defined  by 


D(Ao)  =  ]  e  /  w  e  D(Aj,) 


iV[  Bt  Cj 
0  Ai 


Likewise,  let  Br 


2  be  defined  by 


If  A2  ^  0,  then  define  Zj,  as 

Zj  =  Zj  L2(-r,0;R2)  =  u)  x  r2  x  L2(-r  ,0;  R^) 


with  inner  product 


<z,z>^  =  <z°,z^>0^  +  ';(;,(s)  ,<j;(s)  >  ds 


where  z  =  (w,x,>{i)^  and  z  =  (w.x,^)^  €  Zj,.  respectively. 

As  for  the  Euler-Bernoulli  beam  we  define  A  on  D(A) 


T  w  j  /ye  D(Aj,) 
D(A)  =  •  X  e  Zj,  /  V  6  Hl(-r,0) 

(t,  /  ^(0)  =  X 


2^  by 
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w 

'  AjW  +  Bp  CjX 

A 

X 

- 

A^x  +  A24i(-r) 

0'(*) 

The  Input  operator  B  :  fi  ->■  2^  is  defined  by 


B  - 


0 

Bi 

0 


(84) 


(85) 


The  following  result  follows  from  Theorem  2  and  standard  estimates  from 
the  theory  of  delay  systems. 


THEOREM  5.  The  operator  B  Is  bounded  and  A  generates  a  C^-semlgroup 

At 

e"  on  the  Hilbert  space  2^. 

Again  it  follows  that  a  state  space  model  for  the  Timoshenko  beam  with 
delayed  actuator  dynamics  is  of  the  form 

z(t)  "  Az(t)  +  Bu  (t)  (86) 

c 

where  A  is  defined  by  (83)  -  (84)  and  B  is  defined  by  (85) . 


ISOTROPIC  RECTANGULAR  PLATE 

We  assume  that  the  control  is  applied  such  that  in  equation  (35) 
f(t,x,y)  ■  <|i(x,y)  u^(t)  where  i|>(x,y)  is  a  given  distribution  belonging 
to  H^O)  and  u^(t)  modulates  the  applied  force. 

Let  P  denote  the  "energy  space" 

V  -  H2(n)  n  Hi(n)  (8?) 

0 

where  H^(n)  and  H^(i2)  are  the  standard  Sobolev  spaces  (see  Ref.  16). 

Let  X  ■  L2(0)  and  define  the  differential  operator  K  with  domain 
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D(K)  5  X  by 


0(K)  =  13 

and 

K<^ 


(88) 


(89) 


It  follows  that  the  plate  equation  (35)  with  boundary  condition  (38) 
may  be  written  as  the  second  order  equation 

u^^  +  K^u  *  ij»(x,y)  u^(t)  (90) 

in  L2(0).  We  rewrite  this  as  a  first  order  system  by  using  the  "state" 

■u(t,x) 

z(t)  -  (91) 

_  U^(t,x)  . 


<z,z>  -  D  [Au(x,y)  Au(x,y)]  dxdy 

■'n 

+  m  v(x,y)  v(x,y)  dxdy  (93) 

•'fl 

where  z  ■  (u,v)  and  z  ■  (u,v)  e  2^,  respectively.  As  for  the  beams, 

<z,z>  ■  Hz  11^  is  the  mechanical  energy  at  the  state  z. 

Let  D(A  )  c  2  be  defined  by 
P  -  P 

D(A  )  -  D(K^)  D(K)  (94) 

P 

and 
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A 

P 


■  0  I  ■ 
-K2  0 


(95) 


Also,  define  6 


R 


Zp  by 


B 


0 

L  <i>(x»y) 


(96) 


The  following  Is  a  direct  consequence  of  standard  results  found  In 
Refs.  17,  18  and  19. 


THEOREM  6.  The  operator  B  Is  bounded  and  &  generates  a  C  -group 

P  P  o 

At 

e  on  the  Hilbert  space  2  • 


As  In  the  case  of  the  beam  structures,  we  have  a  state-space  model  for 
the  plate  of  the  form 

z(t)  -  Ap2(t)  +  BpU(t)  (97) 

where  Ap  and  Bp  are  now  defined  by  (94)  -  (96)  . 

GENERAL  FRAMEWORK 

All  of  the  problems  formulated  above  fall  Into  a  specific  class 
of  distributed  parameter  systems  of  the  form 

z(t)  -  A2(t)  +  Bu^(t)  (98) 

where  A  generates  a  C^-semlgroup  on  an  Hilbert  space  2  snd  B  :  R  ->■  Z 
Is  a  bounded  linear  operator.  The  space  Z  Is  called  the  state-space 
and  given  an  Initial  state  zq  e  D(A) >  we  know  that 


Z(t)  »  e^^  Zq  + 


i; 


eA(t-s) 

c 


(99) 
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define  (at  least)  weak  solutions  of  the  partial-differential  equations 
that  govern  the  motion  of  the  particular  structure.  The  particular 
model  (Euler-Bernoulli  beam,  Timoshenko  beam,  Euler-Bernoulli  beam 
with  actuator  dynamics  and  no  delays,  etc.)  determine  the  space  E  and 
the  operators  A  and  B.  We  summarize  these  spaces  and  operators  for  the 
cases  presented  above. 

Euler-Bernoulli:  No  Actuator  Dynamics 
E  =  y  as  defined  by  (41)  -  (42) 

A  =  Ag  as  defined  by  (43)  -  (44) 

B  =  B„  as  defined  by  (45) 
h 

Timoshenko:  No  Actuator  Dynamics 

E  =  ID  as  defined  by  (53)  -  (54) 

A  =  Aj  as  defined  by  (55)  -  (56) 

B  =  B^  as  defined  by  (57) 

Euler-Bernoulli:  Actuator  Dynamics,  No  Delay 
E  =  Ep  as  defined  by  (63)  -  (64) 

A  =  Ag  as  defined  by  (65)  -  (66) 

B  =  Bg  as  defined  by  (67) 

Euler-Bernoulli:  Actuator  Dynamics  with  Delays 
E  *  2g  as  defined  by  (69)  -  (70) 

A  “  A  as  defined  by  (71)  -  (72) 

B  B  as  defined  by  (73) 
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Timoshenko;  Actuator  Dynamics,  No  Dela^ 


2  “  2.^  as  defined  by  (76)  -  (77) 

A  *  Ag  as  defined  by  (78)  -  (79) 

B  “  Bg  as  defined  by  (80) 

Timoshenko;  Actuator  Dynamics  with  Delays 
2  ■  Ej,  as  defined  by  (81)  -  (82) 

A  *  A  as  defined  by  (83)  -  (84) 

B  "  B  as  defined  by  (85) 


Isotropic  Rectanimlar  Plate 


2*2  as  defined  by  (92)  -  (93) 
P 

A  ■  A  as  defined  by  (94)  -  (95) 
P 

B  =  Sp  as  defined  by  (96) 
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THE  CONTROL  PROBLEM 


The  abstract  state-space  formulation  of  the  various  beam  and 
plate  systems  provides  a  convenient  framework  in  which  to  formulate 
and  approximate  control  problems.  Recall  that  in  all  cases  the  number 
of  controls  available  is  finite,  indeed  for  the  problems  considered 
here  we  have  a  single  control. 

As  noted  previously,  the  purpose  of  active  control  is  to  increase 
the  effective  damping;  that  is,  to  enhance  the  decay  rate  of  unwanted 
vibrations.  The  decay  rate  can  be  understood  in  terms  of  the  eigenvalues 
of  the  closed-loop  system.  Over  the  past  thirty  years  a  considerable 
literature  has  developed  on  the  use  of  linear-quadratic  control  theory 
to  design  controllers  for  finite-dimensional  system  (see,  for  example, 
the  excellent  book  by  Kwackernaak  &  Sivan  (Ref.  20)).  While  this  theory 
has  deficiencies,  such  as  requiring  accurate  knowledge  of  the  open- 
loop  dynamics,  it  is  a  useful  part  of  a  comprehensive  control  design 
strategy.  Accordingly,  we  shall  formulate  a  linear-quadratic  control 
problem  for  our  (infinite-dimensional)  Hilbert-space  setting.  The 
control  laws  we  develop  for  the  beam  and  plate  models  are  based  on 
this  formulation. 

All  of  the  state-space  formulations  developed  in  the  section  above 

required  a  Hilbert  space  5,  a  dynamic  operator  A  with  D(A)  c  £  and  an 

input  operator  B:R  -*•  £.  We  assume  that  A  generates  a  C^-semigroup 
At 

e  and  B  is  bounded. 

In  order  to  formulate  an  LQR  problem,  we  introduce  an  output 
operator 
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and  assume  that  C  is  bounded  and  linear.  It  may  be  helpful  to  think 
of  C  as  a  mapping  that  "reads  out"  physical  meaningful  quantities  at 
certain  locations  on  the  structures  (e.g.  velocity  or  strain  measure¬ 
ments).  Mathematically,  we  require  that  C  has  the  form 


(100) 

where  each  "sensor"  i  =  l,2,...,n  is  a  bounded  linear  functional 

from  -  to  .  Examples  will  be  given  below. 

One  can  now  describe  the  control  system  which  is  governed  by  the 
state-equation  in  2 

z(t)  =  Az(t)  +  Bu  (t)  (101) 

c 

with  initial  data 


C  = 
z 


C.z 


G2Z 


C  z 

n 


z(0)  =  ZO  e  d(A)  (102) 

and  output  in  B 

y(t)  *  Cz(t)  (103) 

It  should  be  noted  that  "output"  is  used  here  in  the  sense  that  it  is 
a  variable  to  be  controlled.  The  control  law  we  construct  will 
require  state  feedback  and  we  do  not  consider  the  estimation  problem. 

The  optimal  control  problem  on  J  is  to  choose  a  cc  'trol 
u^  6  L2(0,+“’;B)  to  minimize  the  cost  functional 
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J(2  .U) 


(104) 


I  {ily(t)ll|  +  Rllu^(t)l|2}  dt 

where  S>0isakxk  semi-definite  matrix,  R  >  0  and  llyil^  « 

”  S 

<Sy,y>  is  the  standard  weighted  semi-norm.  It  should  be  noted  that  the 

cost  functional  can  be  generalized  to  include  cross-terms  in  y  and  u 

T 

by  formulating  a  quadratic  weight  on  the  composite  vector  (y,u)  . 

Let  H  :  2  -*■  2  denote  the  self-adjoint  bounded  operator  defined 
by 

Q  -  C*  S  C  (105) 

]|Q 

where  C*  :  8  ->■  2  is  the  adjoint  of  C.  Note  that  one  can  rewrite  the 
cost  function  (104)  in  the  equivalent  form 

1 

J(zo,u)  -  -r  {<(5z(t),z(t)>  +  <Ru  (t),u  (t)>}  dt  (106) 

^  Jo  c  c 

and  hence  the  optimal  control  problem  described  by  (101)  -  (104)  is  the 
direct  generalization  of  the  standard  finite  dimensional  (LQR) -problem. 
Much  of  the  theory  for  the  finite-dimensional  problem  carries  over  to 
this  Hilbert  space  setting.  In  particular,  we  shall  make  use  of  the 
variation  of  parameters  formula 

z(t)  ■  e^*"  zq  +  [  e^^^  Bu  (s)  ds  (107) 

■'0 

for  the  mild  solution  of  (101)  -  (102)  . 

A  function  u  6  L2(0»+";  8)  Is  admissible  control  for  the  initial 
state  £,  or  simply  an  admissible  control  for  if  J(z,u)  is  finite; 
l.e.,  if  the  state  z(t)  corresponding  to  the  control  u(t)  and  the 
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initial  condition  z(0)  =  zq  is  in  L^CO,  <»;£). 

Let  the  operators  A,  B,  Q,  and  R  be  as  defined  above.  A  bounded 
linear  operator  il  :  2  ^  2  is  a  solution  of  the  operator  Riccati 
equation  if  0  maps  the  domain  of  A  into  the  domain  of  A*  and  satisfies 
the  Riccati  equation 

A*n  +  HA  -  n  br'^  b*  n  +  q  =  0  .  (los) 

The  following  result  may  be  found  in  Ref.  15. 

THEOREM  7.  There  exists  a  non-negative  self-adjoint  solution  of  the 

operator  Riccati  equation  (108)  if  and  only  if,  for  each 

z  €  2.  there  is  an  admissible  control  for  the  initial  state 

z.  If  n  is  the  minimal  non-negative  self-adjoint  solution 

of  (108),  then  the  unique  control  u  (•)  which  minimizes  J  and 

c 

the  corresponding  optimal  trajectory  z(*)  are  given  by 

u  (t)  =  -  R~^  B*  n  z(t)  =  -  K  Z(t)  (109) 

c 

and 

z(t)  =  5(t)  z  (110) 

where  S(t)  is  the  C^-semigroup  generated  by 
[A  -  BR~^  B*  n] . 

At 

The  semigroup  e  is  called  the  open-loop  semigroup  and  the  semigroup 
5(t)  is  called  the  optimal  closed-loop  semigroup. 

Note  that  if  the  open-loop  system  is  exponentially  stable,  i.e. 
there  exists  M^,  aj  >0  such  that 

||e^*'||  <  Mj  e"^*^  ,  (111) 
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then  u  (t)  H  0  is  an  admissible  control  for  all  z  e  Z-  In  this  case 
c 

Theorem  7  applies. 

In  the  beam  and  plate  problems  formulated  in  the  previous  sections, 
(111)  does  not  hold.  This  is  due  to  two  factors;  there  is  no  damping 
in  the  models  and  for  the  beam  problems  there  are  rigid  body  modes.  In 
order  to  overcome  these  problems  we  introduce  damping  into  the  models. 

In  particular,  if  equation  (35)  is  replaced  by 

m  u^^(t,x,y)  +  Y  A  u^(t,x,y)  +  D  u(t,x,y)  =  f(t,x,y)  (112) 

where  the  damping  term  y  A  u^(t,x,y)  is  added,  y  >  0,  then  the  basic 
state-space  Zp  =  D  >:  I.2(Q)  is  unchanged.  However,  the  operator  (no 
damping)  Ap  given  by  (95)  does  change  and  for  this  case  (112)  leads  to 
a  new  state  operator 


0  I 

-  K2  -  y  a 


(113) 


-  At 

with  the  same  domain  D(A  )  =  D(A  ).  Moreover,  e  P  will  satisfy  a 

P  P 

decay  rate  of  the  form  (111).  A  similar  (but  slightly  more  complex) 
thing  occurs  if  the  damping  term  y  A  u^(t,x,y)  is  replaced  by  the  so- 
called  Kelvin-Voigt  damping  term  y  A^  u^(t,x,y). 

In  the  beam  problems  we  added  a  viscous  damping  term  to  the  beam 
models.  As  in  the  example  above,  the  state-space  2  does  not  change  but 
the  state  operator  becomes  modified  much  like  Ap  is  changed  to  Ap- 
However,  one  still  does  not  obtain  an  estimate  of  the  form  (111) 
because  of  rigid  body  modes.  It  happens  that  even  though  these  modes 
do  not  have  exponential  decay  rates,  they  are  stabilizable .  Thus,  the 
basic  Theorem  7  can  be  used  to  show  that  for  each  z  6  2  the  LQR  problem 
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for  all  of  the  structures  considered  here  has  a  unique  optimal  solution 
given  by  (109)  where  n  solves  (108) ,  provided  that  we  add  damping  to 
the  structural  elements. 

Another  feature  that  is  apparent  from  the  Riccati  equation  (108) 
is  that  the  adjoint  operator  A*  plays  an  essential  role  in  the  determi¬ 
nation  of  n.  For  the  structural  control  problems,  the  operator  A  is  a 
differential  operator  and  for  the  case  where  there  is  no  damping  and  no 
delays  in  the  actuator,  then  A  *  -  A*.  However,  as  indicated  above  it 
is  essential  to  Include  damping  and  hence  for  all  the  damped  systems  it 
follows  that  A  -  A*. 

The  key  point  is  that  if  one  wishes  to  use  approximations  of  the 
state  operator  A  to  solve  the  Riccati  equation  (108),  then  these  ap¬ 
proximations  must  also  provide  good  estimates  of  A*.  This  is 
particularly  Important  when  delays  are  Included  in  the  model  (or  if  other 
types  of  damping  models  are  used,  see  Refs.  2,  13  and  21).  This  turns 
out  to  be  a  crucial  point  that  is  often  overlooked  in  finite  element 
and  modal  control  approaches  to  these  problems. 

CONSTRUCTION  OF  THE  OUTPUT  OPERATOR 

We  indicate  the  form  of  the  output  operator  for  the  beam  and 
plate  control  systems.  We  provide  details  for  the  Euler-Bernoulli  beam 
problems  and  Indicate  the  necessary  changes  for  the  Timoshenko  models . 

Recall  the  problem  with  no  actuator  dynamics.  For  the  Euler- 
Bernoulli  beam  the  state  is  given  by  (40),  i.e. 
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z(t) 


’  Zl (t) 

0(t) 

Z2(t) 

(jj(t) 

Z3(t) 

n(t)  +  Laj(t) 

Zl4(t) 

C(t)  +  u.(t) 

Z5(t) 

u  (t,x) 

XX 

.*6(t)  . 

_  u^(t,x)  +  xa)(t)  . 

(114) 


Consider  the  problem  where  we  measure  3(t),  ri(t)  +  Lu(t), 

C(t)  +  a)(t) ,  the  average  strain  in  neighborhoods  of  selected  points 
0  ^  Xi  <  *2  ^  *p  1.  1*  ^long  the  beam  and  the  average  velocity  in 


neighborhoods  of  selected  points  0  _<  xj  <  X2  < 
beam. 


< 


X 

q 


£  L  along  the 


We  proceed  now  to  give  the  precise  mathematical  definitions 
needed  to  describe  the  above  output  operator.  Let  e  >  0  and 
X  e  [0,L]  be  a  fixed  point  along  the  beam.  We  consider  the  cases 

e  e 

i)  X  ■  0,  ii)  X  ■  L,  and  ill)  0<x<L.  Ifx-0,  let  m  (x)  »  ]D  (0) 
denote  the  bounded  linear  functional  defined  from  12(0, L)  to  R  by 


[«%)]  4.(0  -  [fll®(0)]  4,(0  f  Ms)  ds  .(115) 

JO 

6  6 

If  X  “  L,  let  n  (x)  *  01  (L)  denote  the  bounded  linear  functional 
defined  by 


[ni®(x)]  <M-)  -  (a*(L)]  «?(•)  -  \  4»(s)  ds 

^L-e 


(116) 


If  0  <  X  <  L  we  pick  e  >  0  such  that  0<x-e<x<x+e<L  and 
let  fll®(x)  denote  the  bounded  linear  functional  defined  by 
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[n^(x)]  ({((O  i,(s)  ds  (117) 

■^x-6 

If  ni*(x)  is  defined  by  (115)  -  (117)  then  the  mean  strain  in  the 

t-neighborhood  of  a  point  x.  with  0<x.  -e<x.  <x.+6<Lis 

till 

given  by 


•'  V  — S 


(118) 


Let  e>0,  0<xi<xo<...<x  <L  and  0  <  Xi  <  Xo  <  ... 

1  _  p  I- 

<  x^  <  L  be  such  that  0  <  x^  -  €  <  x^  <  x^  +  e  <  L,  i  *  l,2,...,p  and 

0  <  X.  -  €  <  X.  <  X.  +  €  <  L,  i  -  l,2,...,q.  Let  C  :  2  ->■  and 
ill  s 

C  :  Z  -►  be  defined  by 

V 

■  [ni*(5i)]  Z5(.)  ■ 

[ffl®(x2))  Z5(.) 

C  z  =  (119) 

9  • 

.  [in®(Xp)]  Z5(-)  J  , 


and 

[*3^(xi)]  zs(’) 

(in®(x2)]  Z6  (• ) 

C  z  =  (120) 

V  • 

_  (x^)  ]  Z6  (• )  .  . 

respectively.  We  now  define  the  output  operator  G„  :  Z  -*■  8  where 

£ 

k»4  +  p  +  qby 
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z. 

i 


^2 

*3 


(121) 


In  summary  then,  for  the  Euler-Bernoulli  beam  without  actuator 
dynamics,  the  output 


y(t)  =  Cg  z(t) 

has  its  first  four  components  given  by  yi(t)  =  9(t),  y2(t)  •  i*>(t) , 

yaCt)  =  [u^(t,L)  +  La)(t)]  and  yt^(t)  =  [u^^(t,L)  +  uj(t)].  The  next  p 

6 

components  are  mean  strains  yj^(t)  =  [ m  (x^)]  i  =  l,2,...,p 

and  the  last  of  components  are  the  mean  velocities  y^(t)  = 

[  Bl  (x^)  ]  [u^(t,*)]  +  x^  w(t)  .  We  can  add  weights  to  the  output  by 
selecting  the  k  x  k  matrix  S  to  be  a  diagonal  matrix  of  the  form 

Sg  =  dia  (c^,  cS,  ...,  c2)  (122) 

where  c^,  i  =  l,2,...,k  are  real  constants. 

The  introduction  of  the  delayed  actuator  dynamics  for  the  Euler- 
Bernoulli  beam  problem  leads  to  additional  states  for  the  actuator 
x(t)  »  col  (x^(t),  X2(t))  and  (in  case  A2  J*  0)  the  history  function 
x(.)  *  col  (xj(.),  X2(«))  6  L2(-r,0;  8^) .  To  account  for  these 
features  we  augment  the  output  operator  defined  by  (121)  and  the 
weighting  matrix  defined  by  (122)  to  include  the  two  actuator  states 
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xi(t)  and  X2(t) .  We  did  not  "sense"  the  history  states  (although  this 
is  clearly  possible)  . 

The  Timoshenko  beam  model  is  treated  in  a  completely  analogous 
fashion.  The  output  operator  for  the  plate  problem  is  also  constructed 

A  A 

to  read  out  the  mean  velocities  at  selected  points  ^ 

i  ■  l,2,...,q  and  mean  values  for  the  Laplaclan  [u^(t ,x,y)  + 

u^(t,x,y)]  at  selected  points  (Xj^>y^)  e  i  ■  1,2 . p.  The  output 

operator  ;  2^  is  again  a  bounded  linear  map . 

FORM  OF  THE  GAIN  OPERATOR  K 


As  noted  above  the  optimal  control  is  of  the  form 

5  (t)  -  -  r"^  B*  n  i(t)  -  -  K  i(t) 
c 

where  K  :  2  -»■  8  is  the  optimal  gain  operator.  The  particular  form  of 
K  depends  only  on  the  choice  of  the  state-space  2  and  the  corresponding 
inner  product  <  ,  >.  In  each  of  the  applications  above,  the  state- 
space  has  a  specific  Hilbert  space  structure  and  the  Rlesz  Representa¬ 
tion  Theorem  (see  Ref.  22)  can  be  applied  to  yield  a  specific  form  for 
the  optimal  feedback  gain  operator.  We  present  these  representations 
for  the  basic  problems. 
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where  k5(*)  e  L2(0,L)  and  eL2(0,L)  are  the  functional  gains  for 

strain  and  velocity,  respectively. 


Timoshenko:  No  Actuator  Dynamics 

Kj.  z^(t)  -  ki  0(t)  +  k2  a)(t) 


+  k3  [u^(t,L)  +  Lw(t)] 

+  +  w(t)] 

rL 

+  kgCx)  [u  (t,x)  +  xuj(t)]dx 
+  [  kg(x)  t<l»^(t,x)  +  uj(t)]dx 

h 

+  [  k7(x)  t4»^(t,x)ldx 

Jo  * 


(124) 


[u^(t,x)  -  \|;(t,x)]dx  . 


Euler-Bemoulll:  Actuator  Dynamics,  No  Delay 

Here  z  ■  z°  ■  col  (z  ,  xj ,  X2)  where  z_  is  defined  by  (40)  and 

Cl  £•  £ 

Z°(t)  -  Kg  Zg(t)  +  k7Xi(t)  +  keX2(t)  .  (125) 

Euler-Bemoulli:  Actuator  Dynamics  with  Delays 

n  T 

Here  z  -  col  (z°,  <)>)  -  (Zg,  xi ,  X2,  ,  <l>2(‘))  and 

,0  fO 

Ked  z(t)  -  KgQ  z°(t)  +  k9(s)  (|>i(s)  ds  +  kio(8)  (>2(3)  ds.  (126) 

•'  -r  ■’  -r 


Timoshenko:  Actuator  Dynamics,  Ho  Delay 

Here  z  ■  z®  ■  col  (z^,  xi ,  X2)  where  z^  is  defined  by  (75)  and 
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Kto  '*■  *^10*2  (t) 


(127) 


Timoshenko:  Actuator  Dynamics  with  Delays 


Here  z  -  col  (zO,  j,)  ■  (z  ,  ,  X2 ,  (()i(0»  i)>2(*))^  and 


,0  ,0 

Kxd  2(t)  »  Kjq  zO(t)  +  k^jCs)  <j>i(s)  ds  +  ki2(s)  <j)2(s)  ds.  (128) 


where  k^(x,y)  ■  A  k  (x,y)  with  ki(*,*)  e  D  and  k2(’.»)  «  L2(i^)  . 

Equations  (123)  -  (129)  Indicate  the  basic  computation  problem 
Is  to  find  numerical  schemes  for  approximating  the  various  gains  and 
functional  gains  that  appear  In  the  representations  (123)  -  (129) . 

This  Is  the  objective  of  the  computational  algorithms  described  In  the 
next  section. 
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NUMERICAL  PROCEDURES 


In  order  to  develop  numerical  schemes  for  solving  the  LQR  problem 
defined  by  (101)  -  (104)  one  must  introduce  approximations  and  analyze 
the  convergence.  In  this  section  we  review  the  basic  ideas  behind  the 
approximation  methods  and  scate  a  general  convergence  result  due  to 
Gibson  (Ref.  14). 

In  order  to  obtain  approximate  solutions  to  the  Riccati  operator 
equation  (108),  it  is  necessary  to  approximate  the  operators  A,  A*,  B, 

B*  and  «).  Since  Q  =  C*  S  C,  this  implies  that  one  must  also  approximate 
C  and  C**  for  the  problems  considered  here,  the  most  difficult  aspect 
of  this  problem  is  the  development  of  convergent  approximation  schemes 

for  A  and  A*«  More  to  the  point,  we  need  approximations  of  the 

At  A*  t 

semigroups  e  and  e 

An  approximating  sequence  for  the  control  problem  defined  by 
equations  (101)  -  (104)  is  a  sequence  of  finite  dimensional  subspaces 
^  c  2,  projections  and  linear  operators 

/  1 


i;*  =  B  *  / 


(130) 


satisfying 


III^  II  I  ,  N  =  1,2,  .  .  ., 
I^z  -»-z  for  all  z  e  2  , 


(131) 
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N 

Bu->-  Bu  > 


for  all  u  e  R 


(132) 


=  [C^l*  S  P^z-*-Q2  »  all  z  e  2  J  • 

At  fA^lt 

The  semigroups  e  ,  e^  must  have  the  property 

(133) 

for  all  z  e  2i  uniformly  for  t  in  compact  intervals  (see  Ref.  14).  The 

construction  of  approximating  sequences  for  the  control  problem  (101)  - 

(104)  is  a  problem  in  numerical  analysis  and  approximation  theory.  The 

basic  idea  is  to  approximate  the  (differential)  operator  by  a  finite- 

dimensional  operator  A  (i.e.  using  finite  elements,  finite  differences, 

N 

modal  truncation,  etc.)  and  then  showing  that  A  converging  to  A 

implies  (133).  This  is  a  nontrivial  problem  in  functional  analysis.  In 

N  N 

fact,  the  basic  question  (when  does  A  P  Az  imply  (133))  is  not  yet 
fully  understood.  A  partial  answer  to  the  convergence  question  is 
provided  by  the  now  famous  Trotter-Kato  Theorem  (see  Ref.  23).  Although 
there  are  a  number  of  extensions  of  this  theorem  (see  Theorem  3.1  in 
Ref.  6),  we  state  a  simple  version  that  is  sufficient  for  the  problems 
considered  here. 


THEOREM  8  (TROTTER-KATO  THEOREM).  Let  A  be  the  generator  of  a  strongly 

At  At  6t  N 

continuous  semigroup  e  satisfying  lie  II  ^  Me*^  .  Assume  that  A  is 


a  sequence  of  operators  generating  strongly  continuous  semigroups 
^A  t  satisfying 
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In  terms  of  numerical  analysis,  condition  i)  is  the  stability  require¬ 
ment  and  condition  ii)  is  the  consistency  requirement.  Thus  convergence 
of  the  approximation  scheme  is  dependent  upon  having  a  consistent  and 
stable  numerical  scheme  that  also  satisfies  the  technical  condition  iii) . 
Any  numerical  scheme  that  does  not  satisfy  conditions  i)  -  iii)  may 
not  produce  an  approximating  sequence  for  the  control  problem  (101)  - 
(104).  Moreover,  even  if  the  approximating  sequence  is  such  that  i)  - 

iii)  are  satisfied,  there  is  no  assurance  that  the  dual  semigroups  will 

f  1  N  At  ^ 

converge  strongly,  i.e.  that  e*^”  P  z  -»•  e  z 

The  following  convergence  results  were  first  established  by 

Gibson  (Refs.  12-14),  We  assume  that  we  have  an  approximating  sequence 

=  (2^,  P^,  A^,  B^,  C^)  that  satisfies  (130)  -  (133)  and  consider  the 

approximating  Riccati  equation 

+  n”  A^  -  n  r"^  [B^]*  n  +  =  o  .  (134) 

N 

The  basic  issue  is  the  convergence  of  IT  to  H . 

N 

THEOREM  9.  Suppose  that  Z  satisfies  (130)  -  (133)  and  for  each  N 

N 

there  is  a  non-negative  self-adjoint  operator  FI  that  solves  (134),  If 
there  exist  constants  M,  M>  0,  0  (independent  of  N)  such  that  the 

closed-loop  semigroup 

N  N  N  -1  N  N 

5(t)  «  exp  ([A”  -  b”  R  ^  [b'"]*  o'"]  t)  (135) 
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satisfies 


ll^^(t)ll  <_  ,  t  >  0  .  (136) 

and 

lln^ll  1  M  ,  (137) 

then  the  operator  Riccati  equation  (108)  has  a  non-negative  self-adjoint 
solution  n  and  for  each  z  e  2 

n  pz  Hz  (138) 

and 

N 

5(t)z^5(t)z  (139) 

uniformly  in  t  >  0. 

The  approximating  gain  operators  are  given  by 

(b"^]*  n”  .  (140) 

The  crucial  issue  is  the  convergence  of  to  the  optimal  gain  operators 
K  defined  by  (109) .  Since  our  structural  control  problems  have  a 
finite  number  of  controls  (i.e.  one  input),  the  following  result 
applies  to  our  problems. 

THEOREM  10.  Assume  all  of  the  conditions  in  Theorem  9  hold.  If  the 
number  of  inputs  is  finite,  then 

^  -  Kii  0  (141) 

as  N 

A  proof  of  this  result  may  be  found  in  Ref.  13  and  in  the  pre¬ 
print  Ref.  14.  In  order  to  use  this  result  one  must  show  that  the 
approximation  scheme  one  proposes  to  use  to  solve  the  control  satisfies 
(130)  -  (133)  and  (136)  -  (137).  Usually  it  is  not  difficult  to 


establish  (131)  and  (132)  and  most  of  the  effort  in  terms  of  numerical 
analysis  comes  from  verifying  conditions  (133)  ,  (136)  and  (137) . 

The  Trotter-Kato  Theorem  is  the  basic  tool  needed  to  establish 
(133)  .  The  work  required  to  establish  (136)  and  (137)  is  essentially 
the  same  as  showing  the  approximating  systems  are  all  uniformly  (with 
respect  to  N)  stabillzable  (see  Ref.  13,  Ref.  14  and  Ref.  24).  This 
can  be  difficult.  However,  for  damped  structures  and  many  of  the 
"standard  schemes,"  (136)  -  (137)  can  often  be  shown  to  hold.  We 
turn  now  to  a  brief  description  of  the  approximation  schemes  for  the 
various  applications. 

While  the  approximation  theory  we  employ  is  rooted  in  the  state- 
space  formulation,  it  is  important  to  employ  'physical'  insights  when 
constructing  specific  approximation  schemes.  When  considering  the 
Euler-Bernoulli  beam  model,  for  example,  it  is  helpful  to  keep  in  mind 
that  the  "infinite"  dimensional  aspect  of  the  state-space  is  engendered 
by  the  fifth  and  sixth  components  of  the  state-vector  and  that  these 
are  each  related  to  the  beam  deflection  [i.e.  zs  ~  u  ,  while  zc  -  u  J. 
Thus,  if  one  thinks  of  approximating  the  deflection  by  a  suitable 
combination  of  shape  function 

u(t,x)  -  z:  a^(t)  h^(x) 

then  the  basis  functions  used  to  approximate  Z5  would  be  the  h|^,  while 
those  used  to  approximate  z  would  be  h^.  In  addition,  the  basis 
functions  should  satisfy  the  essential  boundary  conditions. 


47 


EULER-BERNOULLI  BEAM  APPROXIMATION  SCHEME 

We  consider  a  uniform  grid  with  N  ’panels'  on  the  interval  [0,L] . 

To  represent  deflections  of  the  beam  we  employ  cubic  B-slines,  which 

form  an  (N  +  3)-parameter  family  of  Cl  functions  (see  Refs.  25  and  26). 

N 

The  elements  of  this  family  are  denoted  by  i  ■  -1,0,1, ... ,N+1. 

We  shall  abuse  notations  and  suppress  the  dependence  on  the  mesh 
parameter  N.  The  domain  of  the  {B^}  is  the  interval  [0,L] .  Each  of  the 
can  be  constructed  by  stretching  and  shifting  the  argument  of  a 
fundamental  cubic-spline,  B  :  [8]  ^  R.  The  fundamental  spline  is 
supported  on  the  interval  [-2,2]  (i.e.  it's  identically  zero  outside 
this  interval) .  A  graph  of  the  fundamental  spline  is  shown  in  Figure 
3.  To  construct  an  element  B^(*)  one  employs  the  rule 

B^(x)  =  B(N  .  (x/L)  -  i)  .  (142) 

The  essential  boundary  conditions  for  the  deflections  of  the 
Euler-Bernoulli  beam  require  that  the  shape  functions  and  their  first 
derivatives  should  vanish  at  the  root  (x  *  0) .  Thus,  we  are  led  to  use 
combinations  of  the  B-splines  that  satisfy  these  two  conditions.  This 
leads  to  the  (N  +  1)  parameter  family  of  shape  functions 

hi(x)  =  Bo(x)  -  2  [Bi(x)  +  B_^(x)]  (143) 

h^(x)  =  B^(x)  ,  i  -  2 . N+1  (144) 

The  shape  functions  are  used  to  construct  basis  vectors  in  the 
state-space  i(»  y)  as  defined  by  (41)  -  (42)  .  To  represent  the  "beam- 
velocity"  we  use  the  vectors 

-  [0,  0,  h^(L),  h^(L),  0,  h^(.)]^  (145) 

(i  -  1,2 . N+1)  . 
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Amplitude 
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M 


Note  that  each  e^  has  been  constructed  to  satisfy  the  tip  end-conditions. 
The  'beam-stress'  will  be  presented  by  the  vectors 

®N+l+i  “ 

(i  -  1,2 . N+1). 

In  addition  to  the  beam  deflections  one  must,  of  course,  also 
Include  the  hub  motions  in  the  finite-dimensional  model.  For  the  hub 
velocity  [(i)(t)]  we  employ  the  element 

eo  -  [0,  1,  L,  1,  0,  x]'”^  (147) 

while  for  the  hub  position  the  element  is 

e^^  -  [1,  0,  0,  0,  0.  0]’^  .  (148) 

The  subspace  2**  is  generated  as  the  span  of  the  set  (e^)  and  has 
dimensions  (2N  +4),  i.e. 

2**  ■  span  {e^/i  •  -1,0, . . .  »2N+2}  (149) 

where  the  e^  are  defined  by  (142)  -  (148) . 

Following  the  general  procedure  outlined  above  leads  to  the 
Galerkln  approximation  for  the  system  (101)  -  (103)  restricted  to  the 
subspace  The  form  of  this  model  is 


Gg  a*^(t)  -  a^(t)  +  e”  u(t) 

(150) 

y^(t)  -  c”  a^(t) 

(151) 

- 

N 

The  components  of  a  (i.e.  oq . °2N+2^ 

the  co-ordinates 

for  the  representation  of  the  state  in  terms  of  our  basis  for 
That  is. 
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(152) 


(153) 


where  the  scalars  (i  ■  0,1,..., N+1)  are  defined  by 
f 0  »  1/3  +  II.  +  m  +  I  ]/pl3 
fi  -  [m^h^(L)  L  +  I^h^(L)lpL3 
and  the  (N+1)  x  (N+1)  symmetric  matrix  G  is  given  by 
Gi  =  [m^h^(L)hj(L)  +  I^h^(L)hJ(L) 

+  p  h. (x)  h.(x)  dx]/pL3  . 

h  ^ 

I  is  the  identity  matrix  of  order  (N+1)  and 
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(154) 


where  a  *  (El/ pL**)  is  a  frequency  parameter  and 

G2(l.j)  -  L  f  h"(x)  hy(x)  dx  . 

■’0  J 

N 

Eg  is  a  (2N+4)  component  column  matrix  with  all  zeroes  except  for 
(1/ pL  in  the  second  entry . 

Referring  to  equation  (121)  it  is  seen  that  the  matrix  has 
(4  +  p  -f  q)  rows  and  (2N+4)  columns.  This  is  moat  easily  described  by 
considering  the  individual  rows  (or  blocks  of  them).  The  first  output 
is  the  angular  position  9  and  the  corresponding  row  is  given  by 

C^(l,  •)  »  [1,0, 0,0]  ,  (155) 

where  the  block  0  is  an  (N+1)  vector.  The  second  row  (angular  velocity 
w)  is  given  by 

Cg(2,  •)  ■  [0,1, 0,0]  .  (156) 

The  third  output  is  the  tip-mass  velocity  and  its  row  in  is 

E 

Cg(3,»)  -  [0,L,h®*(L)  ,...,hJ^j^(L)  ,0]  .  (157) 

The  fourth  output  is  the  angular  velocity  of  the  tip— mass  and  is  given 
by 

Gg(^,*)  “  [0,1, h  (L)  , . . .  (L)  ,0]  .  (158) 
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The  next  block  of  p  outputs  correspond  to  measurement  of  the  beam  strain 
at  specified  points  (x^)  on  the  beam.  A  typical  row  in  this  block  is 

Cg  (4  +  i,.)  -  [0.0.0.  ®®(xph"(*). 

...  l®(ij)h;j^^(.)l  (159) 

where  ID  (x^^)  is  the  approximation  to  the  delta  function  (see  equation 
(117)  .  The  final  block  of  q  outputs  describe  measurement  of  the  total 

A 

beam  velocity  at  specified  points  (x^)  on  the  beam.  A  typical  row  in 
this  block  is 

Cg(4  +  p  +  i.  •)  *  [0.  x^.  ID®(x^)  hi(*). 

...  ID®(i^)  0]  .  (160) 

TIMOSHENKO  BEAM  APPROXIMATION  SCHEME 

For  the  Timoshenko  beam  the  state -space  E(»ffl)»R’<BxRxR 
X  L2  ^  1-2  1-2  La  (®®®  (53).  (54)).  The  fifth  co-ordinate  is  related 

to  the  beam's  transverse  velocity  [u^]  while  the  sixth  is  related  to 
the  angular  velocity  [41^].  The  seventh  co-ordinate  is  the  curvature 
(ij;  ]  and  the  eighth  is  the  shear  deformation  [u  -  ip]  .  If  one  thinks 
of  approximating  the  functions  u  and  ip  then  it  is  'natural'  to  coinblne 
the  (N+3)  cubic  B-slines  to  enforce  the  essential  boundary  conditions 
on  u  and  [u^(t.O)  ■  4<^(t.0)  -  0].  This  leads  to  the  (N+2)  parameter 
family  of  shape  functions  that  vanish  at  the  left  end: 

qo(x)  *  ®o(x)  “  4  B_j^(x)  (161) 

qi(x)  ■  4  r  Bi(x)  -  Bo(x)  (162) 

q^(x)  i  B^(x)  .  i  -  2.3,....N+1  (163) 
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The  basis  vectors  for  the  lineal  velocity  are 

=  [0.  0.  q^(L),  0.  q^(- )  .  0,  0,  0]'*^  (164) 

1  •  0,1,..., N+1.  Similarly  the  basis  elements  for  the  rotational 
velocity  are 

®n+2+i  "  ‘*1^*^’ 

1  ■  0,1 . N+1.  To  represent  the  ip  distribution  the  basis  elements 

are 

®2N+4+l  “  (166) 

while  the  basis  elements  for  the  u  distribution  is 

X 

»•  «•  0-  »•  “■  “•  ’;<•>'  • 

As  in  the  Euler-Bernoulli  case,  the  hub  motions  require  the  elements 

e_j^  -  [1  0  0  0  0  0  0  0]^  (168) 

and 

eg  »  [0  1  L  1x1  0  O]'^  .  (169) 

The  subspace  generated  by  the  {e^>  defined  by  (161)  -  (164)  is  again 
denoted  by 

2^  ■  span  {e^  /  1  "  -1,0 . 4N+4}  (170) 

and  is  [4N+61  dimensional. 

Again  one  applies  the  general  theory  for  approximating  the  system 
to  produce  a  model  of  the  form 

gJ  x”(t)  -  hJ  x*’(t)  +  u(t)  (171) 

y^(t)  -  c”  x”(t)  (172) 
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N  N 

In  this  case  the  square  matrices  and  Hij,  are  of  order  (4N+6)  and 

N  N 

is  a  compatible  column  matrix.  The  matrix  has  (4N+6)  columns  and 

(4+p+q)  rows. 

EULER-BERNOULLI  WITH  ACTUATOR  DYNAMICS;  NO  DELAY 

In  this  case  the  state  space  2  (“  2®)  is  the  'Euler-Bemoulli' 
state-space  y  augmented  with  the  two  additional  actuator  states  [see 
(63)  -  (64) ] .  Because  of  the  cascade  structure  of  the  system  (62)  it 
is  relatively  simple  to  extend  the  approximations  used  for  the  Euler- 
Bernoulli  beam  with  no  actuator  (i.e.  (143)  -  (144))  to  produce  the 
approximation  for  2^.  Specifically,  each  of  the  (2N+4)  basis  vectors 
defined  by  (145)  -  (148)  can  be  Injected  into  2g  by  augmenting  two 
zeroes  in  the  places  for  the  actuator  states.  The  additional  basis 
vectors  needed  to  cover  the  actuator  states  are  taken  as 

®2N+5  -  (O'  0)  (173) 

and 

®2N+6  -  (0.  0,  1)  ,  (174) 

where  0  is  the  zero  vector  in  y (defined  by  (41)  -  (42)). 

With  these  basis  vectors  in  hand  it  is  clear  that  the  subspace 
is  (2N+6)  dimensional  and  the  standard  procedure  is  used  to  generate 
the  matrices  needed  to  represent  the  approximating  system.  In  fact  the 
cascade  structure  of  the  system  leads  to  a  simple  block  structure  for 
the  matrices.  More  is  said  about  this  feature  in  the  next  section. 

EULER-BERNOULLI  WITH  DELAYED  ACTUATOR  STATES 

The  inclusion  of  delays  in  the  dynamical  model  for  the  actuator 
leads  to  some  complication  in  the  approximation  procedures.  Because 
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the  actuator  Is  linked  to  the  Euler  beam  in  "cascade"  fashion  the 
state-space  2p  has  the  product  structure  described  above  (see  (69)). 

2g  -  a  X  R2  X  L2(-r,0;82)  . 

Moreover,  the  dynamics  of  the  coupled  system  has  the  cascade  structure 


y(t)  -  Ag  y(t)  +  Bg  Ci  x(t)  (175) 

x(t)  -  Ai  x(t)  +  A2  x(t-r)  +  Bi  u^(t)  (176) 

We  develop  the  approximation  for  this  system  in  two  stages. 
Firstly,  using  the  scheme  discussed  above,  one  approximates  the  operator 
Ag  and  hence  the  Euler-Bernoulli  dynamics.  This  leads  to  the  system 

(150)  which  can  be  put  in  normal  form  by  (effectively)  inverting  the 

w  N  N 

Gram  matrix  Define  the  matrices  A„  and  as 

«  EE 


4  ^  ^4^'"  4 
4  "  4  • 


(177) 

(178) 


This  approximate  Euler-Bernoulli  model  when  coupled  with  the  delayed 
actuator  can  be  written  as  the  delay  differential  system 


z**(t) 


<  4 

0  Ai 
0 


z  (t)  + 


0  0 


0  A2 


z  (t-r) 


u^(t) 

c 


(179) 


Here  z^  has  dimension  (2N+6) . 
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The  system  (179)  is  still  "infinite  dimensional"  because  of  the 
history  term  z  (t-r) .  There  are  many  numerical  schemes  for  constructing 
approximations  for  such  delay-differential  systems  [Refs.  6-8,  13]. 

The  averaging  scheme  [see  Ref.  6]  uses  piecewise  constant  functions  with 
a  uniform  grid  (say  M  sub-intervals)  on  the  interval  [-r.O].  This 
leads  to  an  (approximating)  system  of  ordinary  differential  equations. 


N,M.  . 

V  (t) 


^  gN,M 


(180) 


M  Xf  MM 

where  v  ’  is  a  vector  of  dimension  (M+1) (2N+6) .  The  matrix  A^’  is 
given  by 


(181) 


where  again  I  is  the  identity  of  order  (2N+6) .  The  control  matrix 
(M+l)(2N+6)  X  1  is 
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0 


0 


0 


t 


(182) 


while  the  output  matrix  Is 

pN.M 

N  M 

where  I  Is  the  (2x2)  Identity  and  a  Is  a  scalar.  Clearly  C^*  has 

E 

(h+p+q]  rows;  the  last  two  rows  read  out  the  actuator  'states'  Xi(t) 
and  X2(t) . 

There  are  several  Important  observations  that  can  be  made  about 
the  system  (180)  with  dynamical  matrix  (181)  .  First,  the 

has  a  block  structure  that  Is  Independent  of  the  details  In  the  delay 
differential  equation  (179).  Specifically,  the  matrix  operating  on 
*^(t)  (see  (181))  appears  in  the  upper  left  block  of  while  the 

matrix  operating  on  the  delayed  'state'  z**(t-r)  is  in  the  upper  right 
block.  The  main  diagonal  below  the  (1,1)  block  has  (-M/r)  •  I  In 
each  block,  while  the  subdiagonal  has  blocks  with  (M/r)  •  I. 

The  (1,1)  block.  Itself  has  son«  structure  with  the  Euler  approxi¬ 
mation  [A^]  In  the  upper  left  and  the  (2’<2)  actuator  matrix  in  the 
lower  right.  The  term  in  the  upper  right  depends  on  Cn  and  Hn 
from  the  Euler  approximation  [see  (178)]  and  on  the  actuator  state-to- 
torque  matrix  Ci  [see  (61)]. 


0  al 


0 

0 


(183) 
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The  matrix  in  the  upper  right  block  of  A„’  has  the  (delayed) 

£ 

actuator  matrix  A2  as  its  only  non-zero  entry.  Similarly  the  control 
N  M 

matrix  has  the  (1x2)  actuator  control  matrix  Bi  in  its  only  non¬ 

zero  block. 

N  M 

The  output  matrix  C^’  has  zero  blocks  corresponding  to  the 

’history’  ’states’.  That  is,  the  output  depends  on  the  first  (2N+6) 

element  'block'  in  the  (M+1) (2N+6)  vector  v^*^. 

N  M  N  M  N  M 

In  the  task  of  assembling  the  matrices  A^'  *  ^E* 

thinks  of  the  matrices  C^,  and  as  "data"  (describing  approxi¬ 

mate  Euler  beam  dynamics).  Similarly,  the  matrices  Ai,  A2>  Bi  and  Ci, 
and  the  scalar  r,  are  "data"  for  the  actuator.  This  view  is  also 
exploited  when  the  Timoshenko  beam  is  coupled  to  the  actuator. 

Finally,  since  only  the  last  two  components  of  the  (2N+6)  vector 
N 

z  (t-r)  are  really  needed,  it  is  clear  that  considerable  saving  is 

N  M 

possible.  In  particular  one  can  compress  v  ’  to  length  (2N+6)  +  2M 
without  affecting  the  input-output  behavior  of  the  approximating  system 
[see  Ref.  10].  Indeed,  if  A2  has  rank  (t  e  [0,1,2]),  then 
can  be  compressed  to  length  (2N+6)  +  £M.  Of  course,  with  Z  “  0  the  A2 
matrix  is  null  and  the  system  is  identical  to  the  non-delayed  actuator 


case. 


TIMOSHENKO  BEAM  WITH  ACTUATOR  DYNAMICS 

In  view  of  the  discussion  of  the  structure  of  the  delayed  actuator 
model  for  the  Euler-Bernoulli  beam,  it  is  possible  to  describe  the  case 
for  a  Tlmoehsnko  beam  (with  and  without  delays)  very  simply. 
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N  M 

The  matrices  for  this  case  are  of  the  same  form  as  except 

N 

and  are  replaced  by  the  corresponding  blocks 
from  the  Timoshenko  model.  Specifically 

.  (185) 

=  [gJ]  ^  Cl  .  (186) 

and 

(187) 

The  dimension  of  the  system  is  (4N-i-8)  +  where,  as  before,  ^  is  the 
rank  of  A 2. 


,N.M 


T 

0 


al 


that  the  blocks  A^, 


ISOTROPIC  RECTANGULAR  PLATE 

The  state  space  model  for  this  system  is  given  by  equations  (91)  - 
(96).  In  particular,  if  u(t,x,y)  is  the  transverse  deflection  of  a 
point  on  the  plate  then  our  state  is  identified  with 


z(t) 


u(t,x,y) 

u^(t,x,y) 


(188) 


The  (clamped)  boundary  conditions  require  that  u  vanishes  on  the 
boundary  of  the  rectangular  domain  and  also  that  the  Laplaclan 
(u  +  u  )  vanish  there. 

XX  yy 

Again  one  thinks  of  expressing  the  deflection  in  terms  of  shape 
functions,  say  (p(x,y)  .  Because  the  geometry  of  the  boundary  is 
rectangular,  it  is  convenient  to  express  each  shape  function  as  a 
product  of  one-space  variable  functions.  Specifically,  we  write 
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(189) 


qi(x)  =  Bi(x)  -  ®_j^(x) 

(190) 

q^(x)  =  B^(x)  i  =  2,...,N-2 

(191) 

(192) 

The  spatial  variables  x  and  y  are,  of  course,  independent.  Thus,  when 
constructing  the  B.  for  the  x  variable  the  scale  length  is  L  and  the 

•  X 

grid  parameter  is  [see  equation  (142)]. 

For  the  purposes  of  displaying  the  basis  vectors  it  is  preferable 
to  enumerate  the  approximating  functions  by  a  single  index  (say  k) 
rather  than  a  double  index  (i  and  j) .  Hence  we  define 


Pj^(x,y)  -  qj^(x)  q^(y)  (193) 

where  k  =  (j-1)  •  (N^  -  1)  +  i  and  1  e  [1,2, . . . ,N^-1]  and 

j  e  [1,2,...,N  -1],  There  are  M  =  (N  -1) (N  -1)  independent  functions  p,  , 
y  X  y  K 

and  each  satisfies  the  boundary  conditions. 

The  basis  vectors  used  to  generate  ^  are  the  "displacement" 
elements 
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(194) 


i  «  1,2 . M 

and  the  "velocity"  elements 

lO  ,  p^l  (195) 

i  *  1,2,...,M.  Hence  the  subspace  ^  spanned  by  e^  defined  by  (190)  - 
(195)  has  2(N  -  1)  (N  -1)  elements. 

X  y 

With  the  basis  chosen  the  matrices  needed  to  represent  the 
dynamics  on  the  subspace  ^  can  be  constructed  by  the  standard  Galerkin 
procedure.  Specifically,  one  substitutes  the  approximate  form 


M 

2  (t)  -  I  a^;_(t)  e^ 
k»l 


(196) 


Into  the  system 


z  (t)  *  A  z  (t)  +  0  u(t)  . 
P  P  P  P 


(197) 


The  coefficients  are  then  found  from  the  normal  equations  (see  Refs. 


25,  26)  as 


M  ,M  MM  M 

Cp  a  (t)  -  Hp  a  (t)  +  Ep  u(t) 


(198) 


where  is  the  2M-8quare  (symmetric)  Gram  matrix 


G  (i,j)  H  <e  ,  e,>  , 
P  1  J 


(199) 


l^^d.j)  ■  <A  e^,  e^> 


(200) 


Ep(j>  -  <Bp,ej>  . 


(201) 
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The  inner  product  is  given  by  the  equation  (93) . 

M 

The  output  map  is  found  by  substituting  the  approximate 
form  (196)  into  the  output  equation  (103)  .  has  p  +  q  rows  and 
2M  columns.  The  kth  column  is  given  by 


0“'“'  =  c  e,  . 

P  P  k 


(202) 
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JOINT  MODELING 


A  large  space  station  is  generally  constructed  with  various 
structural  members  joined  together  to  yield  the  desired  configuration. 
The  points  at  which  these  members  are  joined  are  called  the  structural 
joints  and  may  vary  widely  in  type  and  configuration.  The  types  of 
joints  of  interest  here  fall  into  two  categories:  those  whose 
characteristics  are  related  to  friction  properties  or  friction  joints, 
and  those  whose  characteristics  are  related  to  the  properties  of  the 
materials  which  compose  the  joint  or  Integral  joints.  Friction  joints 
are  typically  formed  by  two  friction  surfaces  in  contact  under  the 
action  of  a  constant  clamping  force  such  as  two  beams  joined  with  a 
bolt  or  a  rivet.  Integral  joints  are  formed  by  some  bonding  procedure 
such  as  welding  or  by  shaping  the  connecting  portions  such  as 
inserting  a  rod  in  a  hole  or  some  similar  arrangement.  Clearly  some 
joints  (probably  most)  have  characteristics  associated  with  both  types. 
The  purpose  of  this  section  is  to  examine  various  models  which  can  be 
used  to  describe  the  properties  of  the  joints  with  various  levels  of 
fidelity. 

The  purpose  of  a  joint  is  to  maintain  the  geometric  integrity  of 
the  structure  and  at  the  same  time  transmit  the  loads  from  one  member 
to  another.  Typically  a  fastener  is  used  to  keep  the  lineal  dimensions 
of  the  structure  intact  while  friction  is  used  to  keep  the  rotational 
relation  of  the  structural  members  intact.  A  typical  example  of  this 
type  of  joint  would  be  a  nut,  bolt  and  washer  connecting  two  slender 
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beams.  If  the  bolt  were  loose,  the  joint  would  approximate  a  pure 
pinned  joint,  while  if  tight  it  would  approximate  a  rigid  joint. 
Unfortunately  neither  of  these  two  extremes  are  realizable  and  we  must 
settle  for  something  in  between. 

A  considerable  amount  of  experimental  work  has  been  carried  out 
with  the  object  of  determining  the  properties  of  these  joints, 
particularly  for  the  case  of  structures  which  are  vibrating.  Of  general 
interest  is  the  damping  effect  on  the  motion  while  on  the  other  hand 
particular  interest  is  paid  to  the  mathematical  model  of  the  joint 
which  will  predict  the  detailed  force-displaceinent-time  relationships 
and  also  allow  one  to  predict  the  contribution  to  damping. 

The  bulk  of  investigations  pertain  to  friction  type  joints. 

Beards  and  Woowet  (Ref.  27]  present  experimental  results  indicating  that 
the  amount  of  structural  damping  can  be  improved  by  proper  selection  of 
the  clamping  force  in  friction  type  joints.  They  observed  that  a  low 
value  of  this  force  (maybe  too  low  for  practical  application,  i.e.  the 
structure  loses  some  of  its  integrity)  leads  to  the  optimum  damping. 

They  do  not  concern  themselves  with  a  mathematical  model  of  the  joint. 

An  effort  to  validate  an  early  model  of  a  friction  joint  was  made 
by  Richardson  and  Nolle  [Ref.  28].  Their  work  tests  the  Panovko  [Ref. 
29]  model  as  refined  by  Metherell  and  Oilier  [Ref.  30).  Here  a  normal 
force  P  acts  on  the  cylindrical  friction  joint  consisting  of  two 
cylindrical  members  with  a  friction  interface.  An  application  of  an 
external  moment  causes  friction  shear  stresses  to  be  formed  at  the 
interface  and  a  microscopic  slip  displacement  takes  place  between  the 
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two  members.  This  slip  may  take  place  over  all  or  part  of  the  inter¬ 
face  but  gross  sliding  between  the  two  members  is  excluded  in  the  model. 
The  assumptions  are  that  the  interfacial  slip  occurs  for  all  values  of 
the  load  and  begins  at  the  outer  edge  of  the  contact  circle  and,  as 
the  load  is  increased,  the  slipped  region  extends  radially  inward  as 
needed  to  provide  the  required  load.  The  slipped  region  is  a  circular 
annulus.  The  loading  cycle  is  considered  in  three  parts:  initial 
loading  for  a  previously  unloaded  joint  where  the  moment  goes  from  0 

to  M  ,  unloading  of  the  joint  where  the  moment  goes  from  M  to 
max  max 

M  ,  ,  and  finally  reloading  of  the  joint  where  the  moment  increases  to 
min 

M  again, 
max 

For  the  initial  loading  the  inner  bound  of  the  radius  of  the 
slipped  region  is  determined  from  moment  equilibrium: 


M  = 


rR 

2Ttr^  y  P  dr 
a 


(203) 


where  a  =  inner  radius 

R  =  radius  of  friction  joint 
r  =  radius  to  generic  point 
u  =  coefficient  of  static  friction 
P  =  normal  force  (clamping  force) 

Hence 

a3  =  r3  _  3m/2tt  y  P  (204) 

When  M,  the  moment  becomes  large  enough,  a  0  and  gross  sliding 
occurs.  A  relation  between  the  moment  and  the  joint  displacement  is 
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given  by  (see  Ref.  28) 


<J>  = 


uPR 

GS 


a 

R 


1  ^ 

3' 

[rJ 

(205) 


Here  a  is  dependent  on  the  moment  as  in  equation  (204)  and 

1  1^1. 


GS  ^2^2 


where  the  thicknesses  of  each  part  of  the  cylindrical  joint 

between  body  1  and  body  2,  Gj^,G2  are  the  shear  modulus  of  body  1  and 
body  and  are  properties  of  the  material  adjacent  to  the  friction 
interface. 

During  unloading  the  joint  displacement  is  shown  to  be 


.  - 
‘*’u  GS 


- 

^max 

^max] ^  ,  £ 

fbl 

3' 

-2-3 

R 

+ 

-  +  6 

.  R 

Rj 

■  '  llJ 

. 

(206) 


where  a  *  a  when  M  =  M  <  M  slipping  and  b  is  determined  from 

max  max  gross 


fR 


M  -  M  = 
max 


(2mP)  2iTr-dr  =  (R^  -  b^) 


(207) 


where  M  is  the  value  of  the  moment  during  the  unloading  with  the  lower 

bound  of  M  .  . 

min 

Subsequent  reloading  to  again  provides  the  following  joint 


displacement 


,  .  jiPR 

'’r£  GS 


2-3 


‘max 


^max 


+  6 


■^max 


-  2 


b  )  3 

max 


-I 


+  2 


(208) 


where  C  is  determined  from 
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(209) 


M  -  M  ,  .  (r3  -  C3)  . 

min  J 


The  moment-displacement  characteristics  are  given  in  general  bv 

equations  (206)  and  (208)  as  the  moment  i:  cycled  between  and 

M  .  with  equation  (205)  applying  only  for  the  first  quarter  cycle, 
min 

Since  the  equations  depend  upon  the  fact  that  the  load  is  increasing  or 
decreasing,  the  complete  cycle  forms  a  hysteresis  loop  with  the  energy 
dissipation  given  by 


r^Snax 

°  Jm  ^ 

"min 


3S 


Fl  2  i 

^m] 

^  9  1 

'bml 

1 

1  i  •  / 

L  ' 

R  ’ 

J 

T  4 

R 

J 

(210) 


■'m 


where  —  is  determined  from 
K 


M  .M  . 
max  min  J 


1  - 


IrJ  J  ' 


(211) 


One  of  the  chief  objections  to  the  above  model  is  the  apparent 
difficulty  in  obtaining  the  value  of  u  for  the  static  coefficient  of 
the  joint.  Recall  that  this  is  not  the  value  that  occurs  for  gross 
sliding  but  that  associated  with  a  microscopic  slip  displacement. 

The  friction  joint  discussed  in  the  above  problem  has  a  friction 
interface  between  two  bodies  with  each  body  having  some  thickness  and 
a  torsional  shear  modulus  of  G.  The  load  is  applied  to  the  faces  of 
the  joint  away  from  the  friction  interface.  If  the  thickness  of  these 
two  bodies  approaches  zero,  then  the  restriction  of  no  gross  sliding 
becomes  more  difficult  to  satisfy.  In  the  limit  the  only  motion 
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between  the  two  bodies  would  be  sliding.  The  friction  force  which 
opposes  this  relative  motion  is  what  is  normally  called  coulomb 
friction.  The  model  for  this  type  of  motion  is 


F  - 


if  X  -  y  >  0 
if  X  -  y  <  0 


(212) 


where  x  and  y  are  the  velocities  of  the  two  adjacent  surfaces.  If  the 
motion  is  cyclic,  the  force-displacement  curves  form  a  hysteresis  loop 
and  the  energy  loss  per  cycle  can  be  determined.  Analysis  along  these 
lines  is  presented  in  Refs.  31  and  32. 

The  key  feature  in  reference  31  that  pertains  to  this  work  was 
the  observation  that  under  the  assumption  of  coulomb  friction  several 
types  of  motion  within  the  joint  can  occur  depending  upon  the  proper¬ 
ties  of  the  joint,  e.g.  the  value  of  K  which  is  related  to  the  clamping 
force,  and  the  inertial  properties  of  the  driven  member.  If  K  is  low 
and  the  inertia  high,  then  it  is  likely  that  there  is  always  relative 
motion  between  the  two  surfaces  and  equation  (212)  always  holds.  On 
the  other  hand  if  K  is  large  and  the  inertia  low,  there  will  never  be 
slip  and  the  joint  can  be  treated  as  rigid.  Of  course,  it  is  the 
Interesting  case  which  falls  between  these  two  where  there  can  be  slip 
followed  by  a  non-slip  condition  or  vice-versa  which  is  most  likely  to 
occur.  Reference  31  examines  all  three  possibilities  for  the  case 
where  the  input  displacement  is  sinusoidal.  Results  are  presented  .^n 
the  form  of  phase  plane  diagrams  for  a  joint  with  angular  motion. 
Expressions  for  energy  loss  per  cycle  are  given  noting  that  the  first 
case  of  pure  slip  gives  the  maximum  value  of  loss  and  that  the  value 
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of  K  which  is  best  in  terms  of  energy  loss  is  given  by 


K  =  /2  (213) 

where  A  =  the  amplitude  of  the  motion  (angular),  to  =  the  frequency. 


The  maximum  energy  loss  per  cycle  using  the  K  in  equation  (213)  was 
found  to  be 


E 


4lA^to^ 

TT 


(214) 


where  I  is  the  moment  of  inertia  of  the  damper  or  unit  driven  by  the 
friction. 

All  the  previous  investigations  examined  coulomb  type  friction 
of  some  Interface  between  two  joint  members.  In  addition  the  more 
sophisticated  analyses  Included  consideration  of  the  material  properties 
on  each  side  of  the  interface.  In  all  cases  all  motion  and  structural 
deflections  are  assumed  to  be  in  the  same  plane.  In  References  33 
and  34  the  interaction  of  out-of-plane  deflection  and  the  in-plane 
motion  are  considered.  In  particular  slipping  at  the  Interface  occurs 
in  the  zones  where  the  pressure  (related  to  the  flexural  vibrations) 
is  the  lowest.  Using  these  ideas  these  authors  were  able  to  arrive  at 
an  equivalent  viscous  damping  ratio  for  the  out-of-plane  vibrational 
modes.  A  similar  result  was  obtained  in  a  different  manner  in  Ref. 

35.  Here  a  simpler  model  was  used  for  the  normal  force  variation.  It 
was  assumed  the  normal  force  was  proportional  to  a  normal  displacement 
with  a  constant  K]^  .  This  normal  displacement  is  modeled  as  proportional 
to  the  in-plane  displacement  x,  e.g.  y  ■  cx.  The  joint  friction  force 
in  the  In-plane  direction  is  then  assumed  to  have  the  following  form: 


71 


F  -  -  y  kj  cjx|  sign  (x)  (215) 

Earlier  references  which  support  these  ideas  are  given  in 
references  36  through  40. 

All  the  previous  works  emphasize  the  slipping  of  the  Interface 
between  the  two  joint  members.  If  such  slipping  does  not  occur,  as  in 
the  case  where  the  normal  force  is  extremely  large  and  the  inertial 
loads  are  small,  or  if  the  joints  are  fitted  or  welded  together, 
then  the  joint  properties  become  similar  to  that  of  the  material 
and  can  be  modeled  in  several  different  ways. 

The  simplest  model  is  a  pure  elastic  joint  which  can  be  modelled 
as  a  spring  and  has  no  losses.  Hence  the  force  (or  moment)  trans¬ 
mitted  through  the  joint  is  proportional  to  the  deflection,  F  ■  kx. 

An  equally  simple  model  is  the  viscous  joint  where  the  force  is  pro¬ 
portional  to  the  deflection  rate,  and  not  related  to  the  deflection. 

In  some  cases  slipping  joints  can  be  approximated  by  this  assumption. 

The  force  transmitted  through  the  joint  is  given  by  F  «  cx. 

Two  additional  joint  models  build  on  these  two  by  considering 
the  elements  to  be  in  parallel  (Voigt  model)  or  in  series  (Maxwell 
model) .  The  force  transmitted  would  be  F  =  kx  +  cx  for  the  Voigt 
model.  The  force  transmitted  for  the  Maxwell  model  is  given  by 
F  =  cx^  where 

Xj  -  ^  [x  -  xx]  (216) 

Hence  additional  dynamic  properties  have  been  added  due  to  the  joint 
model . 
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Additional  models  can  be  generated  by  a  combination  of  the  two 
simple  elements  described  above  or  more  conveniently  described  in 
terms  of  the  two  elements,  the  Voigt  and  the  Maxwell  model.  A  four 
element  joint  model  can  be  generated  by  a  Voigt  joint  in  series  with  a 
Maxwell  joint  with  different  spring  and  damping  constants  in  each  part. 
Alternatively  one  could  consider  two  Maxwell  joints  in  parallel,  again 
each  one  with  its  own  spring  and  damping  constants.  Alternative  four 
element  joint  models  are  obtained  by  having  a  spring  element,  a 
damper  element  and  a  Maxwell  element  all  in  parallel,  or  having  two 
Voigt  joints  in  series.  Three  element  joints  can  be  derived  from  the 
four  element  joints  by  letting  one  spring  or  damping  constant  approach 
infinity. 

Typically  these  models  add  dynamics  to  the  joint.  For  a  three 
element  elastic  joint  which  is  characterized  by  a  spring  in  parallel 
with  a  Maxwell  model  (spring  and  damper  in  series) ,  the  transmitted 
force  is  given  by 

F  -  c  Xj  +  kx  (217) 

where 

(x  -  xj)  (218) 

Hence  a  large  variety  of  models  are  possible  to  approximate  a 
joint  made  with  a  viscoelastic  material  and  which  experiences  no 
slipping.  The  models  discussed  can  all  be  represented  by  ordinary 
linear  differential  equations. 

A  more  general  representation  is  suggested  in  Ref.  A1  which  can 
Include  nonlinear  contributions.  For  example  if  a  dead  band  exists 


73 


in  the  joint  a  proposed  force  model  looks  like 


F  -  +  k,x  +  c,;  +  k  x"  +  +  k^lj  +  +  Fj  sign  J 

+  g| xj  sign  (x)  (219) 

where 

■^B  ■  “ 

\b  <*  ■"  -db’ 


and  g  in  the  last  term  corresponds  to  the  appropriate  terms  in  eq.  (215). 
The  generic  form  of  the  force  equation  (219)  is  suggested  to  be 

F  =  k(x,x)x  +  c(x,x)x  +  k  q(c)  (220) 

where 

q(t)  =  Q(x,x)  dt  ,  (221) 

and  Q(*)  represents  some  time  history  of  the  motion  which  contributes 
to  the  "memory"  of  the  material . 

One  form  of  the  memory  term  which  is  linear  that  has  been  suggested 
(see  Ref.  41)  is 

q(t)  =  -  pq  +  x(t)  .  (222) 


The  above  review  represents  the  current  thinking  on  models  to 
represent  the  effects  of  joints  on  the  dynamic  systems.  Most 
researchers  consider  the  models  which  include  local  slipping  the  most 


74 


accurate.  Hence  coulomb  friction  and  its  associated  hysteresis  loops 
would  be  present  in  real  systems.  As  indicated  previously,  these 
effects  have  been  observed  in  the  laboratory  (Ref.  35).  It  would  seem 
appropriate  to  consider  a  general  model  such  as  equation  (220)  and 
fit  the  coefficients  to  match  results  of  laboratory  experiments.  As  a 
last  resort  to  obtain  a  model  which  would  possibly  admit  analysis,  the 
viscoelastic  models  described  by  ordinary  differential  equations  might 
be  the  best  procedure. 

The  inclusion  of  joint  models  in  the  state-space  setting  can  be 
simple  or  difficult,  depending  on  the  choice  of  joint  model.  In  any 
case  numerical  algorithms  have  not  yet  been  developed  for  most  of 
these  models.  The  main  difficulty  with  coulomb  friction  models  is 
that  the  state-space  system  is  a  non-linear  differential  inclusion  that 
requires  new  theoretical  tools  (see  Ref.  19). 
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NUMERICAL  RESULTS 


In  this  section  we  report  on  a  number  of  numerical  experiments 
for  the  control  problems  discussed  in  the  previous  sections.  The  goal 
of  this  section  is  to  demonstrate  the  basic  features  of  the  computa¬ 
tional  algorithms  based  on  the  state-space  models.  The  numerical  runs 
presented  below  are  intended  to  illustrate  the  technique  and  to  point 
out  various  difficulties  associated  with  the  different  models. 

We  conducted  six  runs  on  the  slewing  beam-tip-body  problem  and 
one  run  on  the  plate  control  problem. 

BEAM  PROBLEMS 

For  the  purposes  of  this  report  we  selected  a  beam  (aluminum) 

with  the  following  characteristics; 

L  =  30  ft.  (length  of  the  beam) 

B  ®  1  ft.  (width  of  the  beam) 

h  *  .02083  ft.  (thickness  of  the  beam) 

I^  =  981  slug  ft. ^ 

p  *  5.24  slug/ft • ^ 

m  =  ,327  slug 
c 

I  *  BhVl2  =  7.53520  x  lo"^  ft.** 

E  =  1.44  X  10^  Ib/ft.  ^ 

A  =  Bh  -  .02083  ft.^ 

K'  -  .84746 

G  -  5.76  X  10®  Ib/ft.  2 
I^  *  .3  slug  ft.  ^ 
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Note  that 


El  -  1,085.069  lb  ft.^ 


and 


K'AG  -  10,169,519.98  lb. 


The  control  problem  was  solved  for  the  Euler-Bemoulli  beam  using 
"strain  sensors"  located  at 

Xi“0,  X2“.3L 

and  "velocity  sensors"  located  at 

x^“.5L,  X2=.7L. 

Thus,  for  the  Euler-Bernoulli  beam  there  are  eight  outputs  and  the 
weighting  matrix  is  set  to  be  the  diagonal  matrix 

Sg  *  dia  (100,  100,  102,  102,  102,  102,  100.2,  100.2)  (224) 

For  the  Timoshenko  beam  we  measured  the  mean  values  at  Xj^  ■  0, 

X2  “  .3L  of  27(t)  «  [\)>^(t,x)l  and  Z8(t)  *  [u^(t,x)  -  i);(t,x)l  and  mean 
values  at  Xj  »  .5L,  X2  *  .7L  of  Z5(t)  *  [u^(t,x)  +  xa>(t)]  and  Z6(t)  = 
li|;j.(t,x)  +  oj(t)].  Therefore,  the  weighting  matrix  is  the  diagonal 
matrix 


S^=  dia  (100,  100,  10,  10,  102,  102,  100.2,  100.2,  100,  100,  100,  100)  (225) 

When  actuator  dynamics  were  Included  we  also  "sensed"  the  actuator 
position  x^(t)  and  velocity  X2(t) .  Weights  of  1.0  were  placed  on  these 
states  so  that  for  the  Euler-Bernoulli  beam  with  actuator  dynamics 
the  weighting  matrix  becomes 

Sgjj  -  dia  (Sg,  1.0,  1.0)  ,  (226) 
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and  for  the  Timoshenko  beam  with  actuator  dynamics  the  weighting  matrix 
becomes 

=  dia  (S^,  i.O,  1.0)  .  (227) 

The  weight  on  the  control  is  always  taken  to  be  R  =  1. 

We  considered  actuator  dynamics  with  and  without  delays.  For  the 
model  without  delays  the  system  parameters  were  set  at  (see  equations 
(19)  -  (26)) 

a  I  “  “100  ^12  ~  “  0 

b^  =  10  c^  =  10  Co  =  0 

and  for  the  delayed  actuator  we  used  the  parameters 

a  I  “  “100  ^12  ”  ^21  * 

b^=10  Cj=10  C9  =  0 

Comparison  of  Open  Loop  Systems 

We  started  by  comparing  the  open-loop  systems  corresponding  to 
the  Euler-Bernoulli  and  Timoshenko  beam  models  for  various  lengths  of 
the  beam.  We  used  the  scheme  described  by  equations  (143)  -  (154)  for 
the  Euler-Bernoulli  beam  and  the  scheme  described  by  equations  (161)  - 
(172)  for  the  Timoshenko  beam.  We  denote  by  N  the  number  of  "elements" 
along  the  beam.  More  precisely,  we  sub-divide  the  interval  [0,L]  into 
N  equal  subintervals  and  use  cubic  spline  shape  functions. 

We  added  viscous  damping  to  each  model,  so  that  all  of  the  non¬ 
zero  eigenvalues  had  real  part  equal  to  (-.025).  Both  models  have 
two  zero  eigenvalues  corresponding  to  the  rigid  body  modes.  We  based 


a  22  =  -1.0 
r  =  0.1 


(229) 


a22  =  0 

r  =  0 


(228) 
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calculations  on  N  »  4,  8,  16,  32  and  64  "element"  approximate  models 
and  compared  the  corresponding  open-loop  eigenvalues.  We  let 


x“ 

J 


N 


i 


(230) 


denote  the  jth  eigenvalue  of  the  dynamic  operator  A  . 

In  Table  1  we  show  the  imaginary  parts  of  the  first  eight  eigen¬ 
values  for  an  Euler-Bemoulli  beam  of  length  thirty  feet.  Using  N 
elements  results  in  a  system  of  order  2N  +  4.  Since  there  are  two 
zero  eigenvalues,  we  also  have  N  +  1  complex  conjugate  pairs.  Compar¬ 
ing  the  first  three  columns  of  Table  1  to  the  last,  it  appears  that 
using  N  elements  gives  a  fairly  accurate  approximation  to  the  first 
N  -  2  eigenvalues.  Table  2  shows  the  open-loop  non-zero  eigen- 
frequencies  for  the  Timoshenko  beam  model. 

Observe  that  the  Euler-Bemoulli  model  gives  more  accurate 
estimates  to  the  low  frequency  eigenvalues  and  needs  fewer  elements 
than  the  Timoshenko  model  in  order  to  give  accurate  frequencies.  As 
we  shall  see  later,  this  may  cause  numerical  difficulties  in  the 
control  design  if  the  Timoshenko  model  is  used. 

Recall  that  the  Timoshenko  model  is  more  suitable  for  a  short 
beam  (see  Ref.  3).  We  considered  the  open-loop  eigenvalue  problem 
for  a  1  ft.  beam  to  see  if  there  was  a  significant  difference  between 
the  two  models.  These  results  are  presented  in  Table  3  and  Table  4. 

Tables  3  and  4  show  that  for  the  short  1  foot  beam,  the  Timoshenko 
model  produces  more  accurate  frequencies  than  the  Euler-Bernoulll 
model.  Note  that  the  relative  "error"  between  the  eighth  frequency 
0)0  5  40,244  and  the  16  element  model  prediction  based  on  the  Euler- 
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TABLE  1.  Euler-Bernoulli  Frequencies,  L  *  30  ft. 


N 

4 

8 

16 

32 

64 

N 

U)1 

0.48908 

0.48908 

0.48908 

0.48908 

0.48908 

N 

0)2 

2.1723 

2.1700 

2.1700 

2.1700 

2.1700 

N 

0)3 

6.1752 

6.1077 

6.1046 

6.1044 

6.1044 

N 

0)4 

12.431 

12.077 

12.049 

12.047 

12.047 

N 

27.381 

20.074 

19.419 

19.912 

19.912 

N 

- 

30.129 

29.543 

29.520 

29.519 

N 

0)7 

- 

42.446 

40.714 

40.650 

40.647 

N 

0)0 

- 

57.314 

53.392 

53.238 

53.231 

TABLE 

2.  Timoshenko  Frequencies, 

L  =  30  ft. 

N 

4 

8 

16 

32 

64 

N 

U)1 

0.48961 

0.48935 

0.48914 

0.48896 

0.48884 

N 

0)2 

3.0544 

2.4240 

2.1752 

2.1701 

2.1701 

N 

003 

42.019 

9.7023 

6.249 

6.1065 

6.1045 

N 

(4)4 

278.57 

33.435 

13.250 

12.065 

12.047 

N 

^5 

1219.30 

111.35 

24.925 

19.996 

19.913 

N 

- 

321.49 

43.393 

29.806 

29.522 

N 

0)7 

- 

834.22 

77.314 

41.413 

40.656 

N 

0J8 

- 

- 

144.56 

55.019 

53.252 
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TABLE  3.  Euler-Bernoulli  Frequencies,  L  =  1  ft. 


tl 

4 

8 

16 

32 

64 

N 

0)1 

52.463 

52.463 

52.463 

52.463 

52.463 

N 

0)2 

215.44 

215.44 

215.44 

215.44 

215.44 

N 

0)3 

2,297.20 

2,293.5 

2,293.3 

2,293.3 

2,293.3 

N 

0)4 

6,326.20 

6,218.3 

6,213.8 

6,213.6 

6,213.5 

N 

<^5 

12,496.0 

12,162.0 

12,121.0 

12,119.0 

12,119.0 

N 

- 

20,232.0 

20,001.0 

19,991.0 

19,991.0 

N 

0)7 

- 

30,791.0 

29,867.0 

29,833.0 

29,831.0 

N 

0)8 

- 

44,158.0 

41,746.0 

41,644.0 

41,639.0 

TABLE 

4.  Timoshenko  Frequencies,  L  »  1  ft. 

N 

4 

8 

16 

32 

64 

N 

0)1 

52.462 

52.462 

52.462 

52.462 

52.462 

N 

0)2 

215.31 

215.31 

215.31 

215.31 

215.31 

N 

0)3 

2,343.6 

2,287.6 

2,286.8 

2,286.8 

2,286.8 

N 

0)4 

7,782.1 

6.197.4 

6,173.1 

6,172.8 

6,172.8 

N 

0)5 

28,281.0 

12,260.0 

11,982.0 

11,979.0 

11,979.0 

N 

<^6 

- 

21,440.0 

19,649.0 

19,637.0 

19,637.0 

N 

0)7 

- 

37,089 

29,137.0 

29,083.0 

29,083.0 

N 

0)0 

- 

- 

40,436.0 

40,245.0 

40,244.0 
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Bernoulli  model  oig^  =  41,746  is  3.7%.  On  the  other  hand,  the 
Timoshenko  16  element  model  produces  the  estimate  =  40,436  so 

O 

that  here  the  relative  error  is  less  than  .5%. 

We  turn  now  to  computing  suboptimal  gain  operators.  We  consider 
six  beam  models  and  one  plate  model.  In  each  case  the  goal  is  to 
compute  an  estimate  of  the  optimal  feedback  gain  operator.  We  used 
the  numerical  schemes  and  state  space  formulations  developed  in  the 
previous  sections. 


MODEL  1;  EULER-BERNOULLI  BEAM;  NO  ACTUATOR 
.  State-Space:  2  =  y  defined  by  (41)  -  (42) 

A  =  Ag  defined  by  (43)  -  (44) 
B  *  Bg  defined  by  (45) 

C  =  ^  defined  by  (121) 


Optimal  Feedback  Operator: 

K^Zg(t)  -  Kj  0(t)  +  K2  u)(t) 

+  K  [u  (t,L)  +  Laj(t)] 
3  t 

+  K4  [u^^(t,L)  +  a)(t)  ] 


K5(x)  [u^^(t,x)]  dx 


K6(x)  [u^(t,x)  +  xu(t)]  dx 


(231) 


We  used  the  approximation  scheme  (143)  -  (154)  to  construct  the 
finite  dimensional  system  (150)  -  (151)  and  the  corresponding  approxi¬ 
mate  Riccati  equation  (134).  Potter's  method  (see  Ref.  20)  was  used 
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to  solve  (134)  and  construct  the  approximating  gain  operator 

E 

N 

defined  by  (140).  Again,  K„  will  have  the  form 

E 

N  N  N 

KgZ^(t)  =  Ki  9(t)  +  K2  u.(t) 

+  [Uj.(t,L)  +  Lu)(t)  ] 

+  k'^  +  “(Ol 

N 

+  K5(x)  [u  (t,x)]  dx 

Jq  XX 

N 

+  Ke(x)  (u  (t.x)  +  Xu;(t)]  dx  .  (232) 

•'o 

The  convergence  of  K‘  to  K  is  equivalent  to  the  convergence  of  the 
four  gains 

N  N  N  N 

Ki  -  Ki  ,  K2  -  K.  ,  K3  ^  K3  ,  ^ 

in  R  and  the  two  functional  gains 

K5(x)  -  K5(x)  ,  4(x)  ^  K(3(x) 

in  L2 (0,L) . 

N 

Table  5  shows  the  convergence  of  to  for  i  =  1,  2,  3,  4  and 

N  N 

N  =  8,  12,  16,  32.  Note  that  and  K2  have  "converged  numerically"  for 
N  N 

N  =  16  while  K3  and  K4  seem  to  converge  more  slowly. 


TABLE  5.  Sub-optimal  Gains;  MODEL  1 


N 

4 

4 

8 

-  10.0000 

-  .47911 

4.0716 

-  .24652 

12 

-  10.0000 

-  .47915 

4.0352 

-  .19589 

16 

-  10.0000 

-  .47916 

4.0276 

-  .18769 

32 

-  10.0000 

-  .47916 

4.0155 

-  .18376 
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Figure  4  shows  the  convergence  of  the  functional  gain  KsC*). 

Observe  that  and  are  almost  identical  and  hence  we  do 

N 

not  show  k6‘+(.).  Figure  5  shows  the  convergence  of  KgC*).  Again, 
there  is  no  difference  between  Kf**(*)  and  so  we  plot  only  the 

N  =  32  functional  gain.  Table  6  shows  the  closed-loop  eigenvalues  for 
N  =  8,  16,  and  32.  Recall  that  the  viscous  damping  model  gives  open- 


loop  eigenvalues  X  with  Re(X)  “  -  .025. 


TABLE  6.  Closed-Loop  Eigenvalues;  MODEL  1 


N 

s 

8 

N  =  16 

N  = 

=  32 

-  .046444 

+ 

.047666 

i 

• 

.046444 

+ 

.047666 

i 

-  .046444 

+ 

.047666 

i 

-  .19907 

+ 

.44751 

i 

- 

.19907 

+ 

.44751 

i 

-  .19907 

+ 

.47751 

i 

-  .04155 

2.1694 

i 

- 

.041547 

+ 

2.1693 

i 

-  .041547 

+ 

2.1693 

i 

-  .026758 

+ 

6.1076 

i 

- 

.026758 

+ 

6.1045 

i 

-  .026758 

+ 

6.1043 

i 

-  .025423 

+ 

12.078 

i 

- 

.025419 

+ 

12.049 

i 

-  .025419 

+ 

12.047 

i 

-  .025136 

+ 

20.074 

i 

- 

.025136 

+ 

19.919 

i 

-  .025135 

+ 

19.912 

i 

-  .025109 

+ 

30.129 

1 

- 

.025106 

+ 

29.543 

i 

-  .025106 

+ 

29.520 

i 

-  .025035 

+ 

42.446 

i 

- 

.025029 

+ 

40.714 

i 

-  .025029 

+ 

40.650 

i 

This  example  shows  how  the  computational  scheme  can  produce 
excellent  convergence  for  the  basic  problem  with  an  Euler-Bernoulli 
beam  model.  We  turn  now  to  the  Timoshenko  model. 

MODEL  2:  TIMOSHENKO  BEAM;  NO  ACTUATOR 
•  State-Space:  2  ■  IP  defined  by  (53)  -  (54) 

A  ■  Aj  defined  by  (55)  -  (56' 

B  ■  Bij,  defined  by  (57) 

C  »  Cj.  rfhere 


85 


z  = 


and 


= 


N 

0 

il) 

es 

1 _ 

■  in^(.5L)z5  ■ 

ni^(.3L)z7 

ni^(.7L)z. 

in®(0)zg 

< 

N 

II 

ra®(.5L)z, 

CO 

N 

-3 

-J 

_  I!I®(.7L)zg  _ 

Optimal  Feedback  Operator: 


K^z^(t)  =  Ki  8(t)  +  K2  'j(t) 

+  K3  [u^(t,L)  +  Luj(t)] 
+  Kl,  [  v  (t,L)  +  (>j(t)  ] 


+  I  K^Cx)  (u  (t,x)  +  LoCt'l]  dx 

'0  ^ 


+  1  Kf/x)  (i^  (t,x)  +  M(t)l  dx 
'0 


K^Cx)  [ii;^(t,x)]  dx 


+  1  K 

^0 


g(x)  [u^(t,x)  -  I^t.x)]  dx 


(233) 


(234) 


(235) 
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We  now  use  the  approximation  scheme  defined  by  (161)  -  (172)  to 


construct  the  finite  dimensional  model  (171)  -  (172)  and  the  cor¬ 


responding  Riccati  equation  (134).  Again,  Potter's  method  was  used 


N 

to  compute  the  approximating  gain  operator  and  as  in  the  Euler- 

Bernoulli  case  has  the  same  form  as  K^. 

N 

Table  7  shows  the  convergence  of  to  for  i  =  1,  2,  3,  A. 

N  N  N 

Observe  that  as  in  the  Euler-Bernoulli  beam  Kj,  K2  and  K3  behave 


nicely. 

The  convergence 

of  is  not 

as  rapid. 

TABLE  7. 

Sub-optimal 

Gd^xns  \  MODHLi  2 

N 

K.2 

8 

-  9.9996 

-  .48306 

4.2078 

.55767 

12 

-  9.9996 

-  .48075 

4.2090 

.34931 

16 

-  9.9984 

-  .48025 

4.1900 

9.9691  E-3 

32 

-10.023 

-  .48373 

4.1721 

-8.7446  E-2 

The  convergence  property  of  the  functional  gains  for  the 
Timoshenko  model  is  not  as  well-behaved  as  in  the  Euler-Bernoulli 
model.  At  this  time  we  do  not  have  a  full  understanding  of  this 
problem.  There  may  be  several  reasons  for  this  slow  (or  non-)  con¬ 
vergence,  including  numerical  roundoff  and  the  failure  of  the  Timoshenko 
model  to  approximate  "long  thin"  beams  (see  the  lecture  notes  by 
Taylor,  Ref.  A2) .  Figures  6  through  9  show  plots  for  the  functional 

gains  K^(  •) »  K^(  •) »  K7(*)  and  K0(«).  Note  that  K^(*)  seem  to  remain 

N 

close  for  N  •  8,  12  and  16  and  i  ■  5,  6,  7  and  8.  Only  Kq(»)  seems 
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to  be  converging  in  a  monotone  fashion  for  N  »  8,  12,  16  and  32. 
Although  the  functional  gains  do  not  seem  to  "rapidly  converge",  the 
closed-loop  eigenvalues  do  appear  to  be  correct  in  that  the  N  =  32 
sub-optimal  gain  produces  reasonable  closed-loop  eigenvalues.  This 
is  Illustrated  in  Table  8. 


TABLE  8.  Closed-Loop  Eigenvalues;  MODEL  2 


N  -  8 

N  « 

16 

N 

32 

-.046441 

+ 

.047660 

1  i 

-.046426  + 

.047675 

1  i 

-.046620 

+ 

.047618 

1  i 

-.19916 

+ 

.44774 

i 

-.19934  + 

.44741 

i 

-.19926 

+ 

.44733 

i 

-.042086 

+ 

2.4236 

i 

-.041603  + 

2.1746 

i 

-.042331 

+ 

2.1694 

i 

-.027471 

+ 

9.7022 

i 

-.026820  + 

6.2488 

i 

-.027401 

+ 

6.1064 

i 

-.025743 

+ 

33.435 

i 

-.025451  + 

13.250 

i 

-.025764. 

+ 

12.064 

i 

-.025414 

+ 

111.35 

i 

-.025213  + 

24.925 

i 

-.025232 

19.996 

i 

-.025099 

+ 

321.49 

i 

-.025110  + 

43.393 

i 

-.025061 

+ 

29.806 

i 

-.025069 

834.22 

i 

-.025072  + 

77.314 

i 

-.025067 

+ 

41.413 

i 

MODEL  3:  EULER-BERNOULLI  BEAM;  ACTUATOR,  NO  DELAY 

•  State-Space:  Z  *  2®  defined  by  (63)  -  (64) 

£ 

A  -  Ao  defined  by  (65)  -  (66) 

B  “  Bq  defined  by  (67) 

and 

(236) 


4^E 

*^Eo*E  “  ^Eo 

*1 

- 

*1 

.  *2  . 

.  *2  . 
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Optimal  Feedback  Operator: 


KgZg(t)  +  K^x^(t)  +  KgX2(t) 
9(t)  +  a)(t) 

+  K3  [Uj.(t,L)  +  La)(t)] 

+  K4  [u^^(t,L)  +  U)(t)] 


+ 


fL 

K5(x) 

0 


[Uxj^(t,x)]  dx 


rL 


K6(x)  [u^(t,x)  +  xw(t)]  dx 


+  K  X  (t)  +  K  X  (t)  . 
7  1  8  2 


(237) 


Note  that  the  only  difference  between  MODEL  1  and  MODEL  3  is  the 

addition  of  the  actuator  states  Xj(t)  and  X2(t).  Table  9  shows  the 

N 

convergence  of  the  gains  for  i  »  1,  2,  3,  4>  7  and  8.  Observe  that 
the  convergence  of  each  of  the  finite  sub-optimal  gains  is  much  faster 
than  for  the  case  with  no-actuator  dynamics  (compare  to  Table  5). 


TABLE  9.  Sub-optimal  Gains;  MODEL  3 


N 

k? 

Ks 

8 

-10.000 

-.47902 

4.1250 

-1.4946 

E-2 

-.54392 

-.95773 

12 

-10.000 

-.47902 

4.1241 

-1.5977 

E-2 

-.54392 

-.95773 

16 

-10.000 

-.47902 

4.1242 

-1.6004 

E-2 

-.54392 

-.95773 

32 

-10.000 

-.47902 

4.1241 

-1.6061 

E-2 

-.54392 

-.95773 
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The  addition  of  actuator  dynamics  also  increased  the  rate  of  con- 

N  N 

vergence  for  the  functional  gains  K5(»)  and  K6(*).  This  is  shown  in 
Figures  10  and  11.  In  Figure  11  only  the  graph  of  K^^(*)  is  shown. 

The  graphs  for  N  ■  8,  12,  16  and  32  are  not  distinguishable  on  the 
scale  shown.  The  addition  of  actuator  dynamics  changed  the  gains  as 
expected.  However,  we  did  not  expect  to  see  the  Increase  in  con¬ 
vergence  of  the  functional  gains.  We  do  not  fully  understand  the 
reason  why  the  addition  of  actuator  dynamics  would  improve  the 
numerical  convergence  of  the  optimal  gains. 

As  expected,  the  addition  of  actuator  dynamics  does  affect  the 
location  of  the  closed-loop  eigenvalues.  However,  as  shown  in  Table 
10  the  closed-loop  eigenvalues  of  the  Euler-Bemoulli  beam  with 
actuator  dynamics  are  almost  the  same  as  the  closed-loop  eigenvalues 
for  the  Euler -Bernoulli  beam  without  actuator  dynamics  (see  Table  6) . 
Observe  that  the  second-order  actuator  introduces  two  open-loop  eigen¬ 
values  at 

X,  -  -.50  ±  9.9875  i 
A 

Table  10  gives  the  closed-loop  eigenvalues  for  the  actuator  mode  and 
the  next  seven  frequencies. 
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TABLE  10.  Closed-Loop  Eigenvalues;  MODEL  3 


N  = 

8 

N  = 

16 

N  = 

32 

-5.050 

+ 

8.665  1 

-5.050  + 

8.665  1 

-5.050 

+ 

8.665  i 

-.0463 

+ 

.0475  1 

-.0463  + 

.0475  1 

-.0463 

+ 

.0475  i 

-.1982 

+ 

.4470  1 

-.1982  + 

.4477  i 

-.1982 

+ 

.4477  i 

0 

CM 

0 

1 

+ 

2.169  1 

-.0420  + 

2.169  1 

-.0420 

+ 

2.169  i 

-.0272 

+ 

6.108  1 

-.0272  + 

6.104  1 

-.0272 

+ 

6.104  i 

-.0253 

+ 

12.077  1 

-.0253  + 

12.049  1 

-.0253 

+ 

12.047  1 

-.0250 

+ 

20.074  1 

-.0250  + 

19.919  1 

-.0250 

+ 

19.912  i 

-.0250 

+ 

30.129  i 

-.0250  + 

29.543  1 

-.0250 

+ 

29.520  i 

MODEL  4:  EULER-BERNOULLI  BEAM;  ACTUATOR  DYNAMICS  WITH  DELAY 
.  State-Space:  2  *  2^  defined  by  (69)  -  (70) 

A  =  A  defined  by  (71)  -  (72) 

B  =  B  defined  by  (73) 

C  =  Cgg  defined  by  (236) 

.  Optimal  Feedback  Operator 


fO 

+  K9(s)  4)i(s)ds  + 


:.(s)ds 


Since  we  have  set  d2i  =  0,  one  can  show  that  K 3(3)  =  0  (see  Refs.  6  and 
13).  This  choice  is  made  merely  to  reduce  some  of  the  computational 
burden. 
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In  this  case  the  approximation  scheme  [see  eqns.  (175)  -  (183)] 
Involves  a  spatial  discretization  parameter  [N] ,  and  a  'past-history' 


parameter  [M] .  Based  on  our  earlier  study  of  this  model  [see  Ref.  1], 

we  have  chosen  to  present  results  only  for  the  value  M  *  8.  As  before. 

Potter's  method  was  used  to  compute  a  solution  for  the  Rlccati  equation 

(13A)  and  the  approximating  gain  operator  was  constructed.  The  gain 
N  M 

operator  has  the  form  [see  eqn.  (126)) 


z 


ED 


(t) 


9(t)  +  .(t) 

+  4’“  [Uj.(t,L)  +  Lw(t)] 


+  [  K5***(x)  [u  (t,x)]  dx 

Jq  ** 

+  [  K6’”(x)  [u  (t,x)  +  xa)(t)]  dx 

•'o 

^  „N,M  ,  .  ^  „N,M  ,  . 

+  Ky  Xi(t)  +  Ke  X2(t) 


+ 


(s) 


(<)l(s)  ds 


+ 


rO 


46“(s) 


-r 


(pzCs)  ds 


As  noted  above,  with  a2i  ■  0  in  the  'delay'  matrix  A2  it  can  be  shown 
N  M 

that  Kg*  (s)  *  0.  The  structural  differences  between  this  model  and 
MODEL  3  is  the  functional  gain  Kxo(s)  and  the  dependence  on  two 
discretization  parameters  [N  and  M] . 


lOo 


N  6 

Table  11  shows  the  convergence  of  the  gains  K  ’  for  i  =  1,  2,  3, 

i 

4,  7  and  8.  As  expected,  the  results  are  quite  similar  to  the  MODEL 
3  case  [c.f.  Table  9]. 


TABLE  11.  Sub-optimal  Gains;  MODEL  4 


N 

^N,3 

1 

4'' 

„N,3 

K.3 

„N.8 

„N,a 

8 

-10.000 

-.47888 

4.1044 

-3.3480  E-2 

-.56202 

-.94504 

12 

-10.000 

-.47888 

4.1034 

-3.4480  E-2 

-.56201 

-.94504 

16 

-10.000 

-.47888 

4.1035 

-3.4500  E-2 

-.56201 

-.94504 

32 

-10.000 

-.47888 

4.1034 

-3.4554  E-2 

-.56201 

-.94504 

The  graphs  of  the  functional  gains  K^*^(x)  and  Ke’^Cx)  are  shown  in 
Figures  12  and  13,  respectively.  These  differ  very  little  from  the 
corresponding  results  in  MODEL  3  (see  Figures  10  and  11).  Figure  14  dis¬ 
plays  the  gain  functional  for  the  history  term  K^o"(s).  The  results  for 
N  =  8,  12,  16  and  32  are  indistinguishable.  Note  that  the  history 
variables  have  been  approximated  using  the  "AVE"  scheme  with  M  =  8. 

Thus,  the  graph  shown  in  Figure  14  is  a  piecewise  constant  approxima¬ 
tion  for  Kio(s)  and  the  interval  [-.1,0)  has  been  partitioned  into  eight 
sub intervals . 

The  closed-loop  eigenvalues  for  the  ’actuator  mode’  and  the  first 
seven  'structural  modes'  are  virtually  identical  to  those  shown  in 
Table  10.  It  should  be  noted  that  this  model  also  includes  some 


approximations  for  the  "past  history  modes".  For  this  case  these  have 


significantly  higher  damping  than  the  other  modes. 


MODI’L  5:  TIMOSHKNKO  BEAM;  ACTUATOR  DYNAMICS  WITH  NO  DEI.AY 

•  State-Space:  E  =  s  defined  by  (76)  -  (77) 

A  =  Aq  defined  by  (78)  -  (79) 

B  =  Bg  defined  by  (80) 

C  =  Cg. 

where  Cg  has  the  block  structure 

■  Cp  0  0' 

Cg  =  0  1  0  (238) 

_  0  0  1 

and  Cj.  is  given  in  (233)  -  (234). 

•  Optimal  Feedback  Operator 

+  Kg  [u^(t,L)  +  Lu;(t)] 

+  K4  [.;^^(t,L)  +  w(t)) 

,L 

+  K5(x)  [u^(c,x)  +  xw(t)]  dx 

'0 

+  j  K^(x)  >^(t,x)  +  u)(c)]  dx 

•'0 

-L 

+  K7(x)  [v^(t,x)]  dx 

^0  ^ 

rl 

+  K8(x)  [u  (t,x)  -  ij^(t,x)]  dx 

Jo 

+  Kg  Xi(t)  +  Kio  X2(t)  (239) 
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As  in  the  previous  cases  we  employ  the  approximation  scheme  to 


construct  the  finite-dimensional  model  (185)  -  (187).  Potter's  method 

N 

is  used  to  compute  the  approximate  gain  operator  which  has  the 


same  form  as  [see  eqn.  (239)1 


Table  12  shows  the  convergence  of  the  gain  values  for  i  = 


1.  2,  3,  A,  9  and  10. 

TABLE  12.  Sub-optimal  Gains;  MODEL  5 


N 

K-; 

Kb 

Ki* 

4 

Kio 

8 

-10.000 

-.47759 

3.0496 

-.49331 

-.54523 

-.95785 

12 

-9.998 

-.47927 

4.1594 

.10191 

-.54459 

-.95779 

16 

-9.999 

-.47982 

4.1898 

1.5195  E-2 

-.54488 

-.95782 

32 

-10.002 

-.47992 

4.2168 

1.9758  E-2 

-.54506 

-.95783 

The  first  four  gains  behave  very  much  like  their  counterparts  in  the 
case  without  an  actuator  (see  Table  7].  The  actuator  gains  also 
compare  well  with  the  values  for  the  corresponding  Euler-Bernoulli  case 
[c.f.  gains  and  in  Table  91. 

The  convergence  results  for  the  functional  gains  »  i  ■  5,  6,  7, 
8  are  shown  in  Figures  15  through  18,  respectively.  These  gains  appear 
better  behaved  than  the  case  with  no  actuator  [MODEL  21  but  still  not 
as  good  as  the  Euler-Bernoulli  case  [MODEL  31 .  The  closed-loop  eigen¬ 
values  are  shown  in  Table  13.  The  actuator  mode  is  identical  to 
the  corresponding  Euler-Bernoulli  case  [MODEL  31,  while  the  structural 
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TABLE  13.  Closed-Loop  Eigenvalues;  MODEL  5 


N  = 

8 

N 

= 

16 

N  = 

^  32 

-5.049 

+ 

8.665  i 

-5.050 

+ 

8.665  i 

-5.050 

+ 

8.665  i 

-.0463 

+ 

.0475  i 

-.0463 

+ 

.0475  i 

-.0464 

+ 

.0475  i 

-.1982 

+ 

.4479  i 

-.1985 

+ 

.4477  i 

-.1984 

+ 

.4477  i 

-.0427 

2.4235  i 

-.0421 

+ 

2.1745  i 

-.0420 

2.1694  i 

-.0276 

+ 

9.7022  i 

-.0273 

+ 

6.2488  i 

-.0275 

+ 

6.1063  i 

-.0250 

+ 

33.435  i 

-.0252 

+ 

13.250  i 

-.0253 

+ 

12.065  i 

-.0250 

+ 

111.347  i 

-.0250 

+ 

24.925  i 

-.0250 

+ 

19.996  i 

-.0250 

+ 

321.485  i 

-.0250 

+ 

43.393  i 

-.0250 

+ 

29.806  i 

modes  are  almost  identical  to  the  Timoshenko  beam  without  actuator 
dynamics  (MODEL  2]. 


MODEL  6:  TIMOSHENKO  BEAM;  ACTUATOR  DYNAMICS  WITH  DELAY 
•  State-Space:  c  =  ip  defined  by  (81)  -  (82) 

A  =  A  defined  by  (83)  -  (84) 

B  =  B  defined  by  (85) 

^  ^  ^TD’ 

where 

tpD  =  [Co.  0.  0] 

and  Co  is  given  by  (238). 


•  Optimal  Feedback  Operator 

Kro  ■  Kro 


r 


K^^(s)  (pj(s)  ds  +  ^^(s)  ds  , 


-r 


J-T 
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where  is  defined  by  (239). 

As  In  MODEL  4  we  have  a  cwo-stage  approximaclon  scheme  involving 
both  a  spatial  discretization  parameter  [N],  and  a  past-history 
parameter  [M] .  Again  we  present  numerical  results  for  the  case  M  >  8. 
We  employ  the  approximation  scheme  described  by  equations  (185)  -  (187) 
and  use  Potter's  method  to  solve  the  appropriate  Rlccatl  eqn.  (134). 

The  gain  operator  has  the  form 


0(t)  +  a)(t) 

+  [u^(t,L)  +  La)(t)] 

+  (.|-j.(t,L)  +  u)(t)] 

+  I  (x)  (u  (t.x)  +  xu.(t)]  d: 

+  I  (x)  [<{»  (t,x)  +  a)(t)]  dx 

+  [  Ky’”  (x)  (t,x)]  dx 


J  n 


(x)  [u^(t,x)  -  i;/(t,x)]  dx 


+  K?’”  XI (t)  +  Ki5”  X2(t) 


fO 

-r 

0 


Kli”  (s)  xi(t-s)  ds 


I 

•'-r 


(s)  X2(t-s)  ds 


(240) 


As  was  noted  for  the  case  of  the  Euler-Bemoulli  beam  with  delayed 
actuator  [MODEL  4],  the  special  choice  a2i  ■  0  in  the  delayed  actuator 
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model  implies  that  the  gain  functional  associated  with  the  delayed 

term  xi(t-r)  is  identically  zero.  Thus,  to  approximate  the  gain 

N  M 

operator  we  must  compute  the  six  gain  values  K.’  for  i  =  1,  2,  3,  A, 

N  M 

9  and  10,  and  the  five  functional  gains  K^’  ,  i  =  5,  6,  7,  8  and  12. 

Table  14  shows  convergence  of  the  gain  values  (recall  M  =  8) 
for  N  =  8,  12,  16  and  32.  The  trends  in  this  problem  are  very  similar 
to  the  case  without  the  delay  [c.f.  Table  12]. 


TABLE  14.  Sub-optimal  Gains;  MODEL  6 


N 

4’® 

CO 

00 

^N,8 

K-g 

/.2 

*^10 

8 

-10.000 

-.47771 

3.1485 

-.42911 

-.56338 

-.94517 

12 

-10.001 

-.47926 

4.1686 

.10434 

-.56277 

-.94511 

16 

-10.000 

-.47968 

4.1756 

-3.028  E-4 

-.56300 

-.94513 

32 

-10.003 

-.48239 

4.2755 

8.047  E-2 

-.56669 

-.94547 

The  spatially  distributed  gains  are  shown  in  Figures  19  through 

22.  These  exhibit  the  same  behavior  we  have  seen  in  the 

previous  Timoshenko  beam  models  (MODELS  2  and  5] .  The  history  gain 

Ki2  is  shown  in  Fig.  23  and  is  virtually  identical  to  the  correspond- 

ing  gain  functional  for  the  Euler-Bernoulli  model  with  delayed 
N 

actuator  [see  K[q(s)  in  Fig.  14]. 

While  the  spatial  gains  do  not  behave  very  well  for  this 
Timoshenko  model  the  closed-loop  eigenvalues  seem  to  be  converging 
nicely,  especially  for  the  lower  frequencies. 
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TABLE  15.  Closed-Loop  Eigenvalues;  MODEL  6 


N  = 

8 

N 

a 

16 

N 

= 

32 

-5.526 

+ 

9.064  i 

-5.526 

+ 

9.064  i 

-5.526 

+ 

9.064  i 

-.0463 

+ 

.0475  i 

-.0463 

+ 

.0475  i 

-.0463 

+ 

.0476  i 

-.1982 

+ 

.4480  i 

-.1985 

+ 

.4477  i 

-.1984 

+ 

.4476  i 

-.0425 

+ 

2.4235  i 

-.0419 

+ 

2.1745  i 

-.0427 

+ 

2.1695  i 

-.0275 

+ 

9.7022  i 

-.0272 

+ 

6.2488  i 

-.0280 

+ 

6.1064  i 

-.0250 

+ 

33.435  i 

-.0252 

+ 

13.250  i 

-.0253 

+ 

12.065  i 

-.0250 

+ 

111.347  i 

-.0250  + 

24.925  i 

-.0250 

+ 

19.996  i 

-.0250 

+ 

321.485  1 

-.0205 

+ 

43.393  i 

-.0250 

+ 

29.806  i 

Note  from  Table  15  that  the  actuator  mode  has  somewhat  higher  damping 
and  frequency  than  for  the  Timoshenko  model  without  delay  [MODEL  51; 
they  are  also  higher  than  either  of  the  relevant  Euler-Bernoulli  models 
[MODEL  3  and  4].  The  structural  eigenvalues  are  virtually  identical  to 
the  previous  Timoshenko  modes  [see  Tables  8  and  13]. 

MODEL  7:  ISOTROPIC  PLATE 

The  model,  as  described  previously,  can  be  written  as  the 
second-order  partial  differential  equation 

iii  u^^(t,x,y)  =  D  a‘-  u(t,x,y)  +  u(t)  Kx,y)  .  (241) 

Here  in  is  the  plate  mass,  per  unit  area,  so  that  m  =  P*h  where  p  is 
the  mass  density  of  the  plate  material  and  h  is  the  plate  thickness. 

The  stiffness  parameter  D  is  given  by 

n  r  E  h^  ■ 

°  “  12  (1  - 
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where  E  is  Young's  modulus  for  the  place  material  and  v  is  Poisson's 
ratio.  The  symbol  denotes  the  Laplacian  operator.  The  function 
describes  the  (fixed)  distribution  of  the  applied  control  force  and 
u(c)  is  the  control.  The  plate  has  length  a  in  the  x  direction  and 
b  in  Che  y  direction.  Internal  damping  is  modelled  in  the  Kelvin- 
Voigt  form.  Thus*  the  right-hand  side  of  (241)  is  modified  by  the 
addition  of  a  term 

For  simplicity  the  parameters  taken  for  the  numerical  calculations 
are  D“m*a*b*l.  The  force  distribution  function  is 

<|>(x.y)  -  1  -  [A  (X  -  1/2)  (y  -  1/2)]  2  ,  (242) 

-4 

and  the  damping  parameter  e  is  10 

Following  the  procedures  outlined  in  Section  IV  the  appropriate 
state-space  model  is 

State-Space:  2  ■  Z  defined  by  (92)  -  (93) 

P 

A  -  A  defined  by  (94)  -  (95) 

P 

B  ■  B  defined  by  (96). 

P 

Before  discussing  numerical  results  for  the  control  problem  we  first 

demonstrate  the  model  with  no  input  (open-loop).  The  approximation 

scheme  is  given  by  equations  (196)  -  (200).  In  all  numerical  results 

we  use  equal  numbers  of  divisions  along  the  x  and  y  plate  axes 

(l.e.(  N  ~  N  ■■  N) .  The  dimension  of  the  approximating  subspace  is 
X  y 

2(N  -  1)2.  With  the  parameters  taken  at  unit  value  the  frequency  of 
undamped  motions  of  the  clamped  plate  can  be  shown  to  be  given  by 
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[see  Ref.  5 ,  p.  495] . 


u)  =  (m^  +  n^j 


(243) 


lUtii  =  • 


TABLE  16.  Open-Loop  Frequencies;  MODEL  7 


m 

n 

N  =  4 

00 

K 

"exact" 

1 

1 

19.7428 

19.7394 

19.7392 

1 

2 

49.5106 

49.3556 

49.3480 

2 

2 

79.2358 

78.9707 

78.9568 

1 

3 

101.2630 

98.7953 

98.6960 

2 

3 

130.7811 

128.4034 

128.3049 

3 

3 

181.9428 

177.8226 

177.6529 

Table  16  shows  a  comparison  of  the  exact  and  computed  frequences  for  N 
4  and  8.  It  should  be  noted  that  due  to  the  geometric  symmetry  the 
"off-diagonal"  modes  (m  n)  come  in  pairs.  Only  one  such  mode  is 
shown  in  Table  16,  thus  with  N  =  4  we  have  actually  estimated  nine 
modes . 

The  final  information  needed  to  define  the  control  problem  is 
specification  of  the  cost  function.  The  standard  quadratic  functional 
in  the  form  (104)  is  used  with  control  weight  R  taken  as  unity, 
output  map  C  given  by 


u(t,x,y) 

{n®(.2,  .2)  [4  u(t,x,y)  1 

_  u^(t,x,y) 

_  ni^(.7,  .7)  [u^(t,x,y)  1 

and  diagonal  weights  S  *  diag  [I.,!.]. 
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The  first  controlled  output  is  thus  the  "strain"  (measured  by  the 
Laplacian  of  the  displacement)  at  x  *  y  =  .2,  while  the  second  output 
is  the  velocity  at  the  point  x  =  y  =  .7.  The  feedback  operator  for 
this  model  is  from  (179) 

K  z  (t)  =  [  Ki(x,y)  A  u(t,x,y)  dx  dy 
+  K2(x,y)  u^(t,x,y)  dx  dy 

•'a 

Following  the  same  procedures  as  for  the  previous  models,  we  employ 

the  approximation  scheme  (196)  -  (202)  to  construct  an  approximate 
N  N 

gain  functionals  Ki(x,y)  and  K2(x,y).  Just  as  in  the  open-loop  fre¬ 
quency  calculation  we  have  used  the  same  discretization  parameter  for 
the  two  spatial  variables  (i.e.  =  N) .  Recall  that  the 

approximating  system  has  dimension  2(N  -  1)^. 

We  employ  Potter's  method  to  construct  a  solution  of  the  Riccati 

equation  and  then  compute  the  functional  gains. 

N 

Figures  24-26  present  the  graphs  of  Ki(x,y)  for  N  «  4,  8  and  12. 

Since  these  functionals  depend  on  two  independent  variables  each  N 

value  is  a  separate  plot.  To  provide  some  comparison  we  show  values 

along  the  line  x  »  y  for  K?(x,y),  N  *  4,  8  and  12  in  Fig.  27. 

Figures  28-30  present  the  results  for  the  second  functional  gain 
N 

K2(x,y),  N  =  4,  8  and  12.  Again  the  section  at  x  =  y  is  compared  for 
all  three  N  values  in  Fig.  31. 

The  closed-loop  eigenvalues  are  shown  in  Table  17.  These  all 
appear  to  be  "converging",  except  for  the  damping  in  the  low 
frequency  mode  which  suddenly  jumps  at  N  *  12. 
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Figure  24.  Functional  Gain  ,  N  °  4  ;  Model  7 
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Figure  29 •  Functional  Gain  K2»  N  «  8  ;  Model  7 
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TABLE 

17. 

Closed-Loop 

Eigenvalues ; 

MODEL  7 

N  = 

4 

N  = 

8 

N  = 

12 

-.274  + 

19.746 

i 

-.269 

+ 

19.743 

i 

-1.767 

+ 

19.664 

i 

-.123  + 

49.511 

i 

-.122 

49.356 

i 

-.122 

+ 

49.349 

i 

-.314  + 

79.236 

i 

-.311 

+ 

78.971 

i 

-.312 

+ 

78.959 

i 

-.513  + 

101.263 

i 

-.488 

+ 

98.795 

i 

-.487 

+ 

98.712 

i 

-.855  + 

130.781 

i 

-.824 

+ 

128.403 

i 

-.823 

+ 

128.320 

i 

-1.655  + 

181.943 

i 

-1.581 

+ 

177.823 

i 

-1.579 

+ 

177.678 

i 
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SOFTWARE  ISSUES 

As  part  of  preparing  the  numerical  studies  reported  earlier  it 
was  necessary  to  develop  the  related  research  software.  Since  the 
MODELS  studied  are  interrelated  [i.e.  Euler-Bernoulli  and  Timoshenko 
beam  models;  actuator  models  with  and  without  delays]  some  effort  was 
devoted  to  organizing  the  software  in  an  efficient  way.  The  objective 
of  this  section  is  to  report  on  this  structure  in  the  set  of  software 
procedures.  We  also  assess  what  additions  and  enhancements  would  be 
needed  to  provide  a  reasonable  environment  for  studying  linear 
quadratic  regulator  synthesis  in  a  class  of  models  including  combina¬ 
tions  of  the  beam  and  plate  models  studied  here. 

Starting  from  the  quadratic  regulator  synthesis  procedure,  the 
fundamental  step  is  the  solution  of  a  rather  large  matrix  Riccati 
equation.  As  noted  previously,  the  procedure  used  is  an  implementation 
of  Potter's  method  (see  Ref.  20)*  The  main  calculations  required  in 
this  method  include  finding  the  eigenvalues  and  eigenvectors  of  a 
Hamiltonian  matrix  and  solving  a  linear  system.  Our  implementation 
employed  standard  IMSL  (Version  10)  routines  for  these  calculations. 

In  order  to  use  the  same  procedure  for  all  the  dynamical  systems 
it  was  helpful  to  think  of  the  Potter  procedure  as  an  'operator' 
which  "transforms"  data  files.  This  idea  of  a  'filter'  or  pipeline 
has  been  engendered  by  growing  popularity  of  the  UNIX  operating 
system.  It  is  very  useful  for  the  simplicity  it  produces  when  a 
variety  of  procedures  are  involved. 
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Specifically  in  the  Potter  procedure  we  envision  a  data  set 
consisting  of  the  system  matrices  A,  B  and  C  in  the  now  familiar 
description 

x(t)  ■  A  x(t)  +  B  u(t) 

y(t)  -  C  x(t)  (245) 

The  code  reads  in  the  matrices,  including  the  integer  variables 
describing  the  size  of  each  matrix.  After  the  calculation  is  complete 
the  code  writes  a  file  consisting  of  the  optimal  gain  values  K.  Other 
information.  Including  closed-loop  eigenvalues,  is  also  written. 

Once  the  gain  matrix  is  known  for  a  certain  approximation  it  is 
still  necessary  to  construct  the  components  of  the  gain  operator  (i.e. 
the  gain  values  and  the  functional  gains,  see  eqns.  (123)> 
for  example) . 

To  understand  this  calculation  it  is  helpful  to  recall  that  our 
finite-dimensional  models  in  the  form  [245]  are  coordinate 
representations  of  infinite-dimensional  systems.  That  is,  the  ’state' 
of  the  system  is  given  by 

z”(t)  -  Z  x^(t)  e”  ,  (246) 

N 

where  the  e^  are  the  basis  vectors  chosen  for  a  particular  representa¬ 
tion  [see  eqns.  (145)  -  (148),  for  example]. 

The  matrix  gains  computed  by  the  Potter  procedure  give  the 
optimal  control  as 

u*^*(t)  -  x(t)  ,  (247) 
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where  the  matrix  K  has  one  row  for  each  control  component  (a  scalar 
for  our  examples) .  In  terms  of  the  state-space  one  has 


N* 

u 


(t) 


N  N,  . 

K  z  (t) 


(248) 


and  the  task  is  to  'compute'  the  operator  k  .  In  our  Hilbert  space 
setting  a  bounded  linear  functional  can  be  represented  as  an  inner 
product  so  that  the  control  can  be  expressed  as 

N*  N  N 

u  (t)  -  <K  .  z  (t)>2N  (249) 


where  K  6  Z  . 

Combining  the  representations  given  by  equations  (247)  and  (249) 
one  has 

X  »  <if,  E  x^  ®i^zN  (250) 

Since  (250)  must  hold  for  all  x  and  since  can  be  represented  in  terms 
N 

of  the  basis  {e^}  as 

Kf*  -  I  8^  (251) 

equation  (250)  leads  to 

g”  8  -  [k”i^  (252) 

N 

where  G  is  the  Gram  matrix 

G”(i,j)  -  <e^  ,  ej>  •  (253) 

Once  the  8  coordinates  have  been  determined  one  uses  the  basis 
elements  to  compute  the  'operator'  For  the  results  presented  in  the 

previous  section,  these  calculations  were  done  on  a  PC  using  the  HATLAB 
software.  This  provided  a  convenient  mechanism  for  generating  the 


135 


plots  for  Che  various  functional  gains. 

As  noted  in  the  Numerical  Procedures  section  the  beam  models 
with  actuator  dynamics  can  be  'naturally’  constructed  from  the 
original  models  and  Che  necessary  'data'  describing  the  actuator. 

Thus,  in  its  final  form  our  "beam"  software  set  has  a  single  procedure 
chat  reads  data  files  describing  the  appropriate  beam  model.  This 
data  is  combined  with  the  actuator  description  to  produce  a  system  in 
the  form  (245) .  This  data  is  written  to  a  file  so  that  it  can  be 
read  by  the  Potter  procedure. 

The  software  set  contains  two  codes  for  producing  the  beam 
models.  Each  procedure  reads  the  same  'data'  file  to  describe  the  beam 
properties.  Since  Che  Timoshenko  model  requires  several  parameters  not 
utilized  by  the  simpler  Euler-Bernoulli  model  this  practice  may  seem 
strange.  It  was  adapted  so  Chat  we  could  be  sure  that  the  two  models 
were  given  consistent  data.  This  greatly  facilitated  the  comparison 
reported  in  the  beginning  of  the  Numerical  Results  Section.  The  two 
beam  procedures,  called  Euler  and  Timo,  produce  data  files  containing 
the  A,  B  and  C  descriptions  needed  by  Potter  (or  by  the  procedure  to 
add  actuator  dynamics) .  The  plate  procedure  also  produces  a  file  of 
this  type.  ^ 

Except  for  Che  final  gain  'operator'  calculations,  all  of  Che 
codes  used  in  Che  study  were  written  in  FORTRAN-77.  The  A,  B,  C 
'system'  files  produced  by  the  beam  and  plate  procedures  can  be  quite 
large.  Since  these  are  most  often  not  read  by  the  analyst,  it  was 
decided  to  write  them  as  unformatted  files.  For  debugging  purposes 
it  was  convenient  to  include  a  simple  procedure  to  generate  a  human 
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readable  file  from  this  data. 

The  seven  FORTRAN  procedures  developed  for  the  numerical  experi¬ 
ments  are: 

EULER:  Assembles  (A,B,C)  system  description  for  Euler-Bernoulli  beam 

model 

TIMO:  Assembles  (A,BtC)  system  description  for  Timoshenko  beam  model 

ACT:  Assembles  (A,B,C)  system  description  for  beam  model  with 

actuator  dynamics 

PLATE:  Assembles  (A,B,C)  system  description  for  isotropic  plate  model 

POTTER:  Reads  (A,B,C)  description,  computes  matrix  Rlccatl  solution 
and  feedback  gain  matrix 

GAIN:  Reads  gain  matrix  and  Gram  matrix,  computes  coordinate 

representation  of  gain  operator 

RSYS:  Reads  unformatted  (A,B,C)  system  description  and  writes  out 

formatted  version  of  the  same 

As  a  final  consideration  we  briefly  discuss  the  possibility  of 
generating  a  research  code  for  use  in  studying  LQR  optimal  feedback 
controllers  for  a  class  of  systems  consisting  of  beam  elements  and 
rigid  bodies  coupled  with  joints  and  actuators.  For  definiteness  we 
consider  only  two-dimensional  planar  models. 

Among  the  first  Issues  to  decide  is  which  'elementary'  structures 
are  to  be  allowed  and  then  precisely  how  are  they  to  be  characterized. 
For  discussion  purposes  we  propose  to  include 
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Euler-Bernoulll  beams 


Timoshenko  beams 

Rigid  bodies 

Force/ torque  actuators 

Since  we  are  suggesting  a  linear  dynamical  model,  it  Is  clear  that 
some  features  of  the  motions  are  not  Included. 

Before  beginning  a  discussion  of  some  of  the  technical  Issues  at 

stake.  It  Is  Important  that  we  make  clear  that  we  are  not  contemplating 

N  N  N 

a  computer  code  that  generates  a  f inlte-dlmenslonal  (A  ,B  ,C  ) 
approximating  model  for  such  a  class  of  systems.  The  experience 
reported  In  Ref.  (2)  should  make  clear  that  such  simulation  models  can 
be  treacherous  when  used  In  control  problems.  Rather  we  have  In  mind 
a  computer  aided  procedure  for  assembling  an  infinite-dimensional, 
state-space  model  from  the  "elementary"  models  for  the  constituent 
pieces  of  the  system. 

The  use  of  signal-flow  graphs  or  block  diagrams  to  construct 
models  Including  feed  forward  and  feedback  elements  are  a  weak 
paradigm  for  the  process  we  have  in  mind.  The  parallel  Is  weak  because 
In  the  present  situation  one  has  to  realize  that  not  all  inter¬ 
connections  of  the  elementary  systems  are  permissible.  This  is  because 
the  interconnections  amount  to  imposing  certain  boundary  conditions 
on  the  adjoining  elements.  It  should  be  clear  from  the  state  space 
modelling  olscussed  in  that  section  that  these  boundary  conditions  play 
a  crucial  role  in  the  formulation.  In  our  approach,  the  boundary 
conditions  have  an  effect  on  the  state  space  2  [see  the  plate  model. 
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eqs.  (87)  and  (92)]  and,  more  commonly,  on  the  generator  A  through  its 
domain  [see  the  Euler-Bernoulli  model,  eqn.  (44)]. 

One  aspect  of  the  problem  is  a  system  theoretic  issue:  given  a 
collection  of  dynamical  systems  [z^,  consisting  of  operators  A^, 
and  and  appropriate  spaces  etc.],  under  what  conditions  does 
an  interconnection  produce  a  dynamical  system?  What  is  the  composite 
state-space?  Loosely,  one  expects  that  a  combination  of  two  systems 
Zj  and  Z2  would  result  in  a  product  space  2  =  2i  22-  The  boundary 
conditions,  however  would  tend  to  make  the  correct  space  smaller,  that 
is  the  space  is  some  proper  subset  of  2.  We  are  proposing  that  one 
could  provide  an  automated  procedure  for  constructing  the  composite 
system  Z. 

One  key  aspect  of  such  a  procedure  is  a  convenient  graphics  en¬ 
vironment  for  the  user  interface.  A  combination  of  icons  and  'pop- 
open'  windows  would  provide  cues  for  the  required  data.  A  number  of 
currently  available  system  modeling  packages  provide  such  interfaces. 

A  second  requirement  is  an  underlying  geometry  manager  that  would 
assemble  the  necessary  mappings  from  the  local  geometry  of  the 
individual  elementary  systems  to  the  'world'  coordinates  of  the 
complete  system.  The  final  need  is  the  procedure  that  constructs  the 
composite  state-space  model.  Since  not  all  system  connections  are 
permissible,  this  module  would  need  special  care  to  present  the  user 
with  suitable  prompts  and  diagnostics  should  an  inadmissible  connection 
be  requested. 

In  summary  we  note  that  the  software  we  have  in  mind  does  not 
simply  construct  a  finite  dimensional  approximation  for  the  system 


(that  is  a  slnnilation  model).  Rather  it  is  our  strong  belief  that 
for  control  problems  one  should  retain  the  true  (Infinite-dimensional) 
system  model  as  long  as  possible.  Numerical  approximations  are  to  be 
introduced  only  at  the  last  step,  just  prior  to  computation. 
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CLOSING  REMARKS 


We  have  constructed  state  space  models  and  computational 
algorithms  for  control  of  seven  structural  vibration  problems . 
Convergence  of  feedback  operators  is  dependent  upon  the  numerical 
schemes  satisfying  (131)  -  (132)  and  (133)  -  (137) .  As  indicated 
above  all  of  these  conditions  can  be  rigorously  established  for  MODELS 
1,  3  and  4.  Conditions  (131)  -  (136)  can  be  proven  to  hold  for  all  of 
the  model  problems  considered  in  this  report.  However,  condition 
(137)  is  more  difficult  to  handle.  In  fact,  numerical  evidence  seems 
to  indicate  that  (137)  may  not  hold  for  the  standard  finite  element 
approximation  of  the  Timoshenko  beam  models  (see  Figures  7,  16  and 
20)  and  for  the  finite  element  approximation  of  the  plate  (see  Figures 
24,  25,  26  and  31).  This  is  an  unsolved  problem.  If  one  can  show 
that  (137)  does  not  hold  for  these  schemes,  then  this  would  raise 
several  practical  problems  Involving  the  use  of  finite  element  models 
in  control  design.  This  is  an  area  of  research  that  needs  further 
study. 

The  results  in  this  paper  show  that  practical  software  tools 
for  control  design  can  be  constructed  for  problems  that  Include 
actuator  dynamics  and  time  delays  as  part  of  the  complete  model. 

This  Integrated  approach  has  many  advantages  and  a  few  limitations. 
However,  there  are  several  practical  and  theoretical  questions  that 
must  be  answered  before  a  "computer  aided  package"  could  be  developed 
into  a  useful  design  tool.  These  questions  are  particularly  Important 
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when  dealing  with  Timoshenko  beams  and  plates.  More  fundamental 
research,  numerical  experimentation  and  laboratory  experiments  are 
needed. 
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